Comparison of GRUAN data products for Meisei iMS-100 and Vaisala RS92 radiosondes at Tateno, Japan

A total of 99 dual soundings with Meisei iMS-100 radiosonde and Vaisala RS92 radiosondes were carried out at the Aerological Observatory of the Japan Meteorological Agency, known as Tateno (36.06^∘ N, 140.13^∘ E, 25.2 m; the World Meteorological Organization, WMO, station number 47646), from September 2017 to January 2020. Global Climate Observing System (GCOS) Reference Upper-Air Network (GRUAN) data products (GDPs) from both sets of radiosonde data for 59 flights were subsequently created using a documented processing programme along with the provision of optimal estimates for measurement uncertainty. Differences in radiosonde performance were then quantified using these GDPs. For daytime observations, the iMS-100 temperature is around 0.5 K cooler than RS92-GDP in the stratosphere, with significant differences in the upper troposphere and lower stratosphere in consideration of combined uncertainties. For nighttime observations, the difference is around −0.1 K, and data are mostly in agreement. For relative humidity (RH), iMS-100 is around 1 % RH–2 % RH higher in the troposphere and 1 % RH smaller in the stratosphere than RS92, but both GDPs are in agreement for most of the profile. The mean pressure difference is ≤0.1 hPa, the wind speed difference is from −0.04 to +0.14 m s⁻¹, the wind direction difference is $\le \mathrm{6.4}{}^{\circ }$, and the root mean square vector difference (RMSVD) for wind is ≤1.04 m s⁻¹.

1 Introduction

The Aerological Observatory of the Japan Meteorological Agency (JMA), called Tateno (located at 36.06^∘ N, 140.13^∘ E, 25.2 m above mean sea level), was established in 1920 and has played a leading role in the operation of all JMA radiosonde stations. The Tateno station was chosen as a candidate site for the Global Climate Observation System (GCOS) Reference Upper-Air Network (GRUAN; Seidel et al., 2009; Bodeker et al., 2016) in 2009 and was certified as a GRUAN site in 2018. The Vaisala RS92-SGP radiosonde (referred to here as RS92; Dirksen et al., 2014) was used for routine observations at the Tateno site from December 2009. While RS92-SGP provides data with a high time resolution, i.e. high-altitude resolution, and highly reliable data and has (GRUAN) data products (GDPs), it was necessary to seek an alternative radiosonde model for operational reasons. This is because the payloads often fall within the greater Tokyo metropolitan region (i.e. highly populated area), and therefore, for safety reasons, the use of lighter instrumentation (the weight of RS92 is 290 g) is necessary. In response to this request, the Meisei RS-11G radiosonde (Kizu et al., 2018; Kobayashi et al., 2019; Hoshino et al., 2022a), whose weight is 85 g, was released in 2012 and started being used at Tateno in July 2013. The GDP for RS-11G was developed and certificated in 2019. Following the release of the Meisei iMS-100 radiosonde (Kizu et al., 2018; Hoshino et al., 2022a), smaller ($\mathrm{55}\times \mathrm{53}\times \mathrm{131}$ mm³) and lighter (40 g), model in 2014, it has been used since September 2017. The temperature and humidity sensors of iMS-100 are basically identical to that of RS-11G, but the relative humidity sensor of iMS-100 has a dedicated thermometer, and the global positioning system (GPS) module is different between the two. RS-11G has also been used in operation at Syowa station (Antarctica) since March 2018, while the iMS-100 has been used at Minamitorishima and other JMA stations since 2017 and at stations of other meteorological service providers and numerous research institutes and universities. Syowa and Minamitorishima are currently candidate GRUAN sites.

Figure 1 gives an overview of the GRUAN data stream. The raw data output from the ground system (dc3cb or MWX files for RS92 radiosonde data and JFMT files for iMS-100 radiosonde data) and meta-data files are submitted to the GRUAN Lead Centre (LC) for archiving, and GRUAN LC issues a data processing ticket file, including the data and meta-data file names, IDs for the measurement and the product, the output file name, and other elements for processing. The data processing centre (PC; Lindenberg for RS92 and Tateno for iMS-100 radiosondes) collects raw data, meta-data, and ticket files from GRUAN LC and performs processing to create GDP files in a NetCDF format. The output GDP file are submitted to GRUAN LC. The GDP files are distributed at the National Centers for Environmental Information (NCEI; ftp://ftp.ncdc.noaa.gov/pub/data/gruan/processing/level2/, last access: 12 October 2022) after quality checking.

Figure 1 GRUAN data file stream.

The evaluation of GDP for RS-11G (RS-11G-GDP.1; Kizu et al., 2019), using GDP for RS92 (RS92-GDP.2; Sommer et al., 2012) in dual flight, was discussed in Kobayashi et al. (2019). In Kobayashi et al. (2019), the RS-11G-GDP.1 temperature is around −0.4 K lower than RS92-GDP.2 in the daytime measurement in the stratosphere, while nighttime measurements generally agree well. The RS-11G-GDP.1 relative humidity (RH) was 2 % RH smaller than the RS92-GDP.2 for 90 % RH–100 % RH, and the RS-11G-GDP.1 was around 5 % RH larger than the RS92-GDP.2 at values lower than 50 % RH, so the RH difference exceeds 2 % RH between 500 and 150 hPa in both the daytime and nighttime data. The pressure difference was 0.5 hPa in the troposphere, and the geopotential height difference was around 10–20 m in the stratosphere.

Although the sensors of iMS-100 are almost identical with RS-11G, the data processing algorithms have some distinct updates for IMS-100-GDP, such as the correction of the hysteresis effect of the RH sensor and the correction of geoid model used in GPS module.

In this paper, Sect. 2 describes the instrumentation used and GRUAN data products (including a brief description of the data processing) for iMS-100 and RS92. Section 3 outlines the dual soundings, and Sect. 4 describes the comparison analysis method. The comparison results are given in Sect. 5, and Sect. 6 discusses outcomes from a dual flight of iMS-100 and a chilled-mirror hygrometer, Meisei SKYDEW (Sugidachi, 2019). Section 7 summarises the findings. The related abbreviations are listed in Appendix A.

2 Instrumentation

2.1 Sensor material and specifications

Table 1 shows the specifications of iMS-100 and RS92. More detailed specifications for each model are described in Vaisala Oyj. (2013), Dirksen et al. (2014), Meisei Electric Co., Ltd. (2020), and Hoshino et al. (2022a). The ground system for iMS-100 is Meisei MGPS2 and that for RS92 was Vaisala DigiCORA III and was replaced with Vaisala MW41 in September 2019.

(Hyland and Wexler, 1983)

Table 1 Specifications of the radiosondes (Vaisala Oyj., 2013; Meisei Electric Co., Ltd., 2020; Dirksen et al., 2014; Hoshino et al., 2022a). Note: GBAS is ground-based augmentation system, SBAS is Satellite-Based Augmentation System, and DWH is depth, width, and height.

^a Wind also can be derived from Doppler shift in GPS signal, like RS-11G. ^b Weight for biodegradable materia model, while 38 g for EPS model.

2.2 GRUAN data processing for iMS-100

Figure 2 gives an overview of GRUAN data processing for iMS-100. The ticket file (^*.gpt) contains the raw data file name (^*JFMT.DAT) and meta-data file name (^*.gmd), the individual identification number for the observation in the GRUAN data archive, and the file name for processed data. GMDB is the abbreviation for the GRUAN Meta-data Data Base, which contains information on payload equipment (such as radiosondes, balloons, parachutes, unwinders, and rigs). In pre-processing, the lag of the data acquisition timing in the processor between temperature/humidity and GPS-related data (time and positioning) are adjusted, relative humidity is recalibrated with additional ground check data (at 0 % RH and 100 % RH conditions, where available), geometric altitude is corrected with the finer geoid model (see Sect. 2.2.3), and the initial data for ascending are extracted by identifying the start and end of ascent. Usually, the end of ascent is when the radiosonde reached its maximum altitude (thus, the start of the descent or the end of radiosonde signals), but ascent data are truncated to ensure their reliability in the event of missing temperature values (i.e. gaps in observation) over a period of 180 s or more. Initial data are processed at each step, as outlined below, to produce the final data. Derived data related to humidity (such as partial water vapour pressure, frost point temperature, and precipitable water vapour) are also calculated. To support climate record quality, GDP data include their uncertainty at each measurement point. The processed data are stored in NetCDF format, and the GDP for iMS-100 radiosonde data is denoted as iMS-100-GDP. The processing algorithm is documented in Hoshino et al. (2022a) for detail, so the outline and difference from RS-11G-GDP.1 (Kizu et al., 2018; Kobayashi et al., 2019) is described in this paper.

Figure 2 GRUAN data processing for iMS-100 (GMDB from GRUAN Meta-data Data Base). The circle is the process starting point. Parallelograms and rectangles with double lines on the left and right sides represent the data and sub-processes for calculation, respectively.

2.2.1 Temperature

An overview of the temperature calculation is given in Fig. 3. T₀, U₀, P_surf, and H_fin are uncorrected temperature values, uncorrected relative humidity values, surface pressure, and corrected geopotential height (see Sect. 2.2.3). If the length of the string or the unwinder is equal to or less than 10 m, then temperature spike filtering is applied when pressure is below 30 hPa. For this filtering and radiation correction based on the heat balance equation (JMA, 1995), provisional pressure P₀ is derived using T₀, U₀, surface pressure data, P_surf, and H_fin (Sect. 2.2.3). T₁, the temperature after spike filtering, is T₀ if no filtering is applied. The ascent rate, asc, at the previous sampling time is used as the ventilation speed, vent, for the calculation of the radiation correction. The ventilation speed at the release time is set as 0 m s⁻¹. The moving averages of T₁ within ±1 s are calculated to determine the uncorrected temperature data with radiation correction, T₂. The amount of radiation correction, T_{rad_corr}, is calculated using T₂, P₀, P_surf, vent, and the solar elevation angle, which is derived from GPS time, ϕ₀ (latitude), and λ₀ (longitude). T_{rad_corr} is subtracted from T₂, to determine T₃. As the thermistor is not completely spherical, there is a small dependence on the angle to sunlight. Accordingly, the Kaiser filter is applied to T₃ to minimise this orientation effect on obtaining the final temperature value, T_fin. The uncertainty budget for temperature is shown in Table 2. The uncertainties are classified as being either correlated or uncorrelated when the source of the uncertainty is systematic or random, respectively (Dirksen et al., 2014). For example, the uncertainties due to the smoothing processes are classified as uncorrelated uncertainties, and uncertainties due to calibration or solar radiation correction are considered to be correlated uncertainties. Figure 4a and b show the vertical profiles of temperature and its uncertainty at 00:00 UTC on 15 September 2017. The dashed line in Fig. 4a indicates the tropopause height. The blue, green, and red lines in Fig. 4b are correlated, uncorrelated, and total uncertainty, respectively. The blue and red lines almost overlap, but the correlated uncertainty is dominant. For daytime observation, the correlated uncertainty increases with altitude because the amount of the solar radiation correction increases, and the uncertainty from it also increases.

Figure 3 Temperature processing for iMS-100 (L is the total string, or unwinder, length; SEA is the Sun elevation angle derived from time and the radiosonde position). The meanings of the shapes are same as in Fig. 2. Diamonds represent the decision.

Table 2 Sources contributing to iMS-100 temperature measurement uncertainty.

Figure 4 The profile of (a) temperature, (b) temperature uncertainty, (c) RH, (d) RH uncertainty, and (e) pressure uncertainty at 00:00 UTC on 15 September 2017. The dashed line in panel (a) is the tropopause height. The blue, green, and red lines in panels (b), (d), and (e) are correlated, uncorrelated, and total uncertainty, respectively. Black lines in panels (b), (d), and (e) show the criteria for screening described in Sect. 4.1.

2.2.2 Relative humidity (RH)

An overview of the RH calculation is given in Fig. 5. T_hum and L are the temperature of the humidity sensor and length of unwinder or string, respectively. For the time lag correction, which corrects the effect of the time constant delay, the response time coefficients are constant for the RS-11G-GDP.1 but are dependent on whether the RH sensor is absorbing or desorbing water vapour molecules (i.e. whether RH is increasing or decreasing, respectively) for the iMS-100-GDP. The trend of RH change is estimated from the slope of the tangent of U₀, using the method proposed by Savitzky and Golay (1964), using the U₀ values within ±5 s values. The response time, τ_hum, is calculated from the obtained trend value together with the temperature of the RH sensor, T_hum. Time lag correction is applied to U₀ to determine U₁ values, which are separated into low- and high-frequency components, U_1low and U_1high, via the digital filtering with a cutoff period set as 4 times the period of pendulum motion, T_pend. High-frequency components sometimes exhibit spike-like variations, which are considered to be associated with water droplets or ice attached on the humidity sensor. These spike-like variations are removed by using a minimum-pass filter, moving averages, and an infinite impulse response (IIR) filter to determine U_2high. The low- and high-frequency components, U_1low and U_2high, are summed to create U₂. For the RS-11G-GDP.1, the hysteresis errors, i.e. the small error remaining after the time lag correction, possibly due to the inhomogeneity of change in the thin-film polymer RH sensor (termed the “slow-regime” in Dupont et al., 2020), are not corrected but only considered to be an uncertainty component. For the iMS-100-GDP, the estimation of the hysteresis errors are formulated as follows: the slow part ratio and its response times for absorption and for desorption are assumed as 3 % and 300 and 12 000 s, respectively. The RH of the slow part, U_delay, and the hysteresis-corrected RH, U₃ are derived via iterative calculation with the following:

$$\begin{array}{}\text{(1)}& {\displaystyle }{U}_{\mathrm{delay}}\left(t\right)={\displaystyle \frac{U_{\mathrm{2}}\left(t\right)+\left(\mathrm{1}-\mathit{\alpha }\right){\mathit{\tau }}_{\mathrm{hys}}{U}_{\mathrm{delay}}\left(t-\mathrm{1}\right)}{\left(\mathrm{1}-\mathit{\alpha }\right){\mathit{\tau }}_{\mathrm{hys}}+\mathrm{1}}}\text{(2)}& {\displaystyle }{U}_{\mathrm{3}}\left(t\right)=\left(\mathrm{1}+{\mathit{\tau }}_{\mathrm{hys}}\right)U_{\mathrm{delay}}\left(t\right)-{\mathit{\tau }}_{\mathrm{hys}}{U}_{\mathrm{delay}}\left(t-\mathrm{1}\right),\end{array}$$

where α is the slow part ratio, τ_hys is the response time for the slow part, and U₃ is the corrected RH.

Figure 5 As per Fig. 3 but for RH.

The temperature–humidity dependence (TUD) correction is applied to minimise RH biases at lower temperatures in the same way as for RS-11G-GDP.1 by Eqs.(3) and (4), but the coefficients are redetermined via additional chamber experiments and flight comparison results with the Cryogenic Frostpoint Hygrometer (Vömel et al., 2007, 2016; Table 3).

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{\Delta }U& =\left(K_{\mathrm{0},\mathrm{0}}+K_{\mathrm{0},\mathrm{1}}\times U_{\mathrm{3}}+K_{\mathrm{0},\mathrm{2}}\times U_{\mathrm{3}}^{\mathrm{2}}\right)\left(\mathrm{1}\right.\\ & \left.-\left(K_{\mathrm{1},\mathrm{0}}+K_{\mathrm{1},\mathrm{1}}\times U_{\mathrm{3}}\right)\right)\exp \left(-T_{\mathrm{hum}}/\left(K_{\mathrm{2},\mathrm{0}}\right.\right.\\ & \left.\left.+K_{\mathrm{2},\mathrm{1}}\times U_{\mathrm{3}}\right)\right)+\left(K_{\mathrm{3},\mathrm{0}}+K_{\mathrm{3},\mathrm{1}}\times U_{\mathrm{3}}+K_{\mathrm{3},\mathrm{2}}\times U_{\mathrm{3}}^{\mathrm{2}}\right),\end{array}\end{array}$$

$$\begin{array}{}\text{(4)}& U_{\mathrm{4}}=U_{\mathrm{3}}-\mathrm{\Delta }U.\end{array}$$

Table 3 Coefficients for TUD correction (Eq. 3).

During soundings, the temperature of the RH sensor is not necessarily equal to the ambient air temperature because of solar heating and the thermal lag of the RH sensor. Values monitored with an RH sensor that are higher (lower) than those of the ambient air would result in dry (wet) biases. Accordingly, U₄ values need to be translated to RH with respect to the saturation vapour pressure associated with air temperature using the following equation (referred as the $T_{\mathrm{s}}/T_{\mathrm{a}}$ correction).

$$\begin{array}{}\text{(5)}& U_{\mathrm{fin}}=U_{\mathrm{4}}\times {\displaystyle \frac{e_{\mathrm{s}}\left(T_{\mathrm{hum}}\right)}{e_{\mathrm{s}}\left(T_{\mathrm{a}}\right)}},\end{array}$$

where e_s(T) is saturated water vapour pressure over liquid water, based on the formulation by Hyland and Wexler (1983), as follows:

$$\begin{array}{}\text{(6)}& \log e_{\mathrm{s}}={\displaystyle \frac{a_{\mathrm{0}}}{T}}+a_{\mathrm{1}}+a_{\mathrm{2}}T+a_{\mathrm{3}}{T}^{\mathrm{2}}+a_{\mathrm{4}}{T}^{\mathrm{3}}+a_{\mathrm{5}}{T}^{\mathrm{4}}+a_{\mathrm{6}}\log T,\end{array}$$

where T is the temperature (K) with the coefficients, a_i, as shown in Table 4.

Table 4 Coefficients for saturated water vapour in Eq. (6).

U_fin is the final value of RH. The uncertainty budget for RH is shown in Table 5. The uncertainties due to the calibration and time lag correction are classified as correlated uncertainties. The uncertainty due to the low-pass filter is uncorrelated uncertainty. The uncertainty due to the IIR filter has correlated components propagating from the time lag correction and uncorrelated components from the low-pass filter. The uncertainties from the hysteresis correction, the TUD correction, and the $T_{\mathrm{s}}/T_{\mathrm{a}}$ correction are calculated by propagation from the uncertainties from all previous steps. The final correlated uncertainty is obtained from synthesising the calibration uncertainty and the correlated uncertainty after the $T_{\mathrm{s}}/T_{\mathrm{a}}$ correction is applied. Figure 4c and d show the vertical profiles of RH and its uncertainty at 00:00 UTC on 15 September 2017. The blue, green, and red lines in Fig. 4d are correlated, uncorrelated, and total uncertainty, respectively. The green and red lines almost overlap at most altitudes. This means that the uncorrelated uncertainty is dominant. The correlated uncertainty has some magnitude (1 % RH–4 % RH) in the troposphere but is nearly negligible above 50 hPa level. If U_fin+u(U) is less than 0 % RH, then these are treated as missing values (NaN) in output.

Table 5 Sources contributing to iMS-100 relative humidity measurement uncertainty.

2.2.3 Geopotential height

An overview of the geopotential height calculation flow is given in Fig. 6. The raw geometric altitude, Z₀, is calculated from ellipsoidal height derived with GPS signal and the internal geoid model of iMS-100, but the grid resolution of the internal geoid model with $\mathrm{10}{}^{\circ }\times \mathrm{10}{}^{\circ }$ is too coarse, and the interpolated geoid height may differ from the actual geoid height by 20 m or more in some regions (near Japan is the one of the regions with large differences; see Fig. 7). Accordingly, geometric altitude is recalculated using the geoid model (Pavlis et al., 2012) with a finer grid (${\mathrm{5}}^{\prime }\times {\mathrm{5}}^{\prime }$) to derive Z₁. The results of the verification using the Global Navigation Satellite System (GNSS) simulator show that Z₀ was found to have a ≃1 s delay with respect to the assumed altitude. Z₁ is corrected for this delay to obtain Z₂. Also, this verification shows that there is a delay of several seconds in measurements just after the launch. This delay is attributed to Kalman filtering in positioning by the GPS module, and the difference becomes negligible several seconds after launch. Accordingly, the altitude in this transition period is interpolated with the known release altitude and the observed altitude at the end of transition to obtain Z₃, which is then applied with a moving average to determine the final geometric altitude, Z_fin. The geopotential height, H_fin, is calculated using the following:

$$\begin{array}{}\text{(7)}& H_{\mathrm{fin}}={\displaystyle \frac{g_{\mathit{\varphi }}}{g_{\mathrm{0}}}}{\displaystyle \frac{R\times Z_{\mathrm{fin}}}{R+Z_{\mathrm{fin}}}},\end{array}$$

where R is the radius of the Earth (6 378 136.0 m), g₀ is the standard gravity acceleration (9.80665 m s⁻²), and g_ϕ is the gravity acceleration at the observation latitude.

Figure 6 As per Fig. 3 but for geometric altitude, Z, and geopotential height, H.

Figure 7 (a) Geoid height for iMS-100-GDP (EGM2008; ${\mathrm{5}}^{\prime }\times {\mathrm{5}}^{\prime }$), (b) iMS-100 original ($\mathrm{10}{}^{\circ }\times \mathrm{10}{}^{\circ }$), and (c) the difference (northwestern Pacific basin).

The geopotential height uncertainty consists of the components due to vertical positioning precision (correlated) and smoothing (uncorrelated).

2.2.4 Pressure

The iMS-100 radiosonde is not equipped with a pressure sensor. In recent years, radiosondes without a pressure sensor have been widely used. This is in part because of reduction in the weight and manufacturing costs. Also, the accuracy of the atmospheric pressure derived using geopotential height obtained from GPS measurements and hydrostatic equation-based approximation is sufficient in the troposphere and superior in the stratosphere (e.g. Nash et al., 2011).

Pressure is calculated from the geopotential height, H_fin (Sect. 2.2.3), temperature, T_fin, and RH, U_fin, using the hydrostatic equation.

$$\begin{array}{}\text{(8)}& P_i=P_{i-\mathrm{1}}\exp \left[-{\displaystyle \frac{g_{\mathrm{0}}}{T_{\mathrm{v}m}{R}_{\mathrm{d}}}}\mathrm{\Delta }H\right],\end{array}$$

where P_i is the calculated pressure (hPa), P_i−1 is the pressure at the previous level (or time step, i.e. 1 s before), ΔH is the thickness (m) between the two consecutive levels, $\mathrm{\Delta }H=H_i-H_{i-\mathrm{1}}$, g₀ is the standard gravity acceleration, T_vm is the average virtual temperature between the two consecutive levels (K), and R_d is the specific gas constant of dry air. The virtual temperature is calculated as follows:

$$\begin{array}{}\text{(9)}& T_{\mathrm{v}}=T{\displaystyle \frac{(\mathrm{1}+r_{\mathrm{v}}/\mathit{\epsilon })}{\mathrm{1}+r_{\mathrm{v}}}},\end{array}$$

with

$$\begin{array}{}\text{(10)}& r_{\mathrm{v}}={\displaystyle \frac{\mathit{\epsilon }e}{P-e}},\end{array}$$

where T_v is the virtual temperature (K), r_v is the mixing ratio of water vapour, ε is the ratio of the molecular weight between water vapour and dry air (=0.622), P is pressure, and e is water vapour pressure calculated from temperature and RH based on the water vapour pressure equation by Hyland and Wexler (1983). For P, P_i−1 is used. The pressure uncertainty is calculated at each data point based on the error propagation rule from temperature, RH, and geopotential height uncertainties from the previous time step. The propagation for correlated and uncorrelated uncertainties is calculated separately and combined to derive the total uncertainty. An example of the pressure uncertainty profile is shown in Fig. 4e.

2.2.5 Wind

An overview of the wind calculation flow is given in Fig. 8. While the initial (unsmoothed) wind speed and direction are derived from GPS Doppler shift for RS-11G-GDP.1, the initial wind speed wspeed₀ and direction wdir₀ are derived as motion vectors from longitude (λ; lon₀ in Fig. 8) and latitude (ϕ; lat₀ in Fig. 8) based on GPS positioning for iMS-100-GDP. The great circular distance, d, and direction, θ, between the position at t_i (λ_i, ϕ_i) and t_i+1 (λ_i+1, ϕ_i+1) are given with spherical trigonometry by the following:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}d& =R\arccos \left(\sin {\mathit{\varphi }}_i\sin {\mathit{\varphi }}_{i+\mathrm{1}}\right.\\ & \left.+\cos {\mathit{\lambda }}_i\cos {\mathit{\lambda }}_{i+\mathrm{1}}\cos \left({\mathit{\lambda }}_{i+\mathrm{1}}-{\mathit{\lambda }}_i\right)\right),\end{array}\end{array}$$

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\mathit{\theta }& =\mathrm{90}\\ & -\arctan \left({\displaystyle \frac{\cos {\mathit{\lambda }}_i\tan {\mathit{\lambda }}_{i+\mathrm{1}}-\sin {\mathit{\varphi }}_i\cos \left({\mathit{\lambda }}_{i+\mathrm{1}}-{\mathit{\lambda }}_i\right)}{\sin \left({\mathit{\lambda }}_{i+\mathrm{1}}-{\mathit{\lambda }}_i\right)}}\right),\end{array}\end{array}$$

where R is the radius of the Earth, and d and θ are rendered as wspeed₀ and wdir₀, respectively. The smoothing process is same as RS-11G-GDP.1. The wind vectors with a wind speed of wspeed₀ and a wind direction of wdir₀ are decomposed to the zonal and meridional wind components, u₀ and v₀. Each of these components are smoothed by a low-pass filter using a Kaiser window with a cutoff frequency of $f_{\mathrm{cutoff}}=\mathrm{1}/T_{\mathrm{pend}}$, for which T_pend is the cycle of the pendulum motion by string with the length L to derive the final wind components of u_fin and v_fin, respectively. The final wind components are synthesised to wind speed, wspeed_fin, and direction, wdir_fin.

Figure 8 As per Fig. 3 but for wind.

Uncertainties in the zonal and meridional wind components consist of those related to the positioning uncertainty and those from smoothing by the Kaiser filter. Both components are uncorrelated uncertainties. The positioning component is estimated as 2.5 m s⁻¹, which is statistically derived using the GPS simulator. The uncertainties in wind speed and direction are calculated from the partial differentiation of the formula that converts wind components to wind speed or direction.

2.3 GRUAN data processing for RS92

As data processing for RS92-GDP version 2 has been detailed by Dirksen et al. (2014), only a brief description is provided here.

Raw temperature data are corrected for solar radiation and heat spike errors. Solar radiation errors relate to the overall direct and scattered solar irradiance, ambient pressure, and ventilation and are estimated at the GRUAN LC from a radiative transfer model that takes into account the solar elevation angle at the time of monitoring. Vaisala radiation error correction data are also available in table form. GRUAN data processing for RS92 involves an application of the average of the two, as it remains unclear which of the two correction models is more appropriate. Heat spike errors are removed via a low-pass digital filter with a cutoff frequency of 0.1 Hz.

RS92 RH sensors have a temperature-dependent dry bias. GRUAN data processing corrects for this based on multiplication with an empirical correction factor before other forms of the correction are applied. Raw RH data are corrected for radiation dry bias, sensor time lag, and temperature-dependence errors. Radiation dry bias is caused by solar heating on the RH sensors, and the same approach as for the temperature sensor is used to estimate the amount of the correction required. The RH sensor response slows at low temperatures, and time lag becomes significant below −40 ^∘C. This is corrected based on the relationship between a time constant and temperature using a low-pass filter in the GRUAN data product for RS92 (Dirksen et al., 2014).

The RS92 used at Tateno has a pressure sensor and a GPS receiver, both of which can be used to calculate geopotential height. Pressure measurement data are used to derive geopotential height in the lower part of the profile where the signal-to-noise performance of the pressure sensor is sufficiently good, and measurements from the GPS sensor are used in the upper part of the profile. The altitude of the switch is typically between 9 and 17 km (Sommer et al., 2016). The pressure sensor is recalibrated against the reference value from a station barometer during the ground check, and a calculation is performed to determine the correction factor for application to the entire pressure profile during sounding (Dirksen et al., 2014). U and V data are retrieved from the Doppler shift in the GPS carrier signal, and noise is removed using a low-pass digital filter. The smoothed data are converted into wind speed and direction values (Dirksen et al., 2014). Uncertainties in each parameter for RS92 are described in Dirksen et al. (2014).

While the authors used version 2 of the RS92 GDP, version 3 is supposed to be available in the near future (Sommer et al., 2016), and it would be useful to redo the analysis with it.

3 Method used for dual sounding

Dual soundings with iMS-100 and RS92 for intercomparison have been conducted from September 2017 to January 2020 at 00:00 UTC (09:00 LT; daytime) or 12:00 UTC (21:00 LT; nighttime) once a week, except for wind conditions in which a payload may fall to populated areas around metropolitan Tokyo (in general, from July to mid-September, with relatively weak westerly winds in the upper troposphere due to the northward displacement of the subtropical jet stream and the stratospheric easterly winds). There were 99 flights during this period (52 daytime and 47 nighttime).

Payload configuration for dual sounding is shown in Fig. 9. A 1200 g balloon was used for all dual soundings. The iMS-100 and RS92 radiosondes were hung on both ends of a 0.9–1 m corrugated plastic board or bamboo rod, and this rig was covered with aluminium tape to reduce the effects of radiation, based on the proposal of von Rohden et al. (2016).

Figure 9 Payload configurations for dual soundings.

Figure 10 The profile of (a) temperature, (b) RH, and (c) the uncertainties for RH at 00:00 UTC on 9 November 2018. The horizontal dashed line in panel (a) represents the tropopause. The blue and red lines in panel (c) show the criteria of screening, according to Eqs. (14) and (16). Panels (d) and (e) show whether the data within (blue) or exceeding (red) the criteria by Eqs. (14) and (16) at each altitude, respectively. Panels (f) and (g) show the percentage of data within or exceeding in the whole profile and are, thus, equivalent to the rearrangement of panels (d) and (e).

4 Method for comparison

GDP data both for iMS-100 and RS92 are collected at 1 s intervals. Temporally simultaneous measurements were compared using the two statistical approaches adopted by Kobayashi et al. (2019) to evaluate differences in the data products. The first method, which is the layer mean in the ensemble of dual flights (described in Sect. 4.4), is to obtain a profile of the differences between the two data products. This is necessary to create a homogenised data set for the climatological discussion at a site where there were instrument changes. The second method, which is the consistency verification using the uncertainty per single measurement level following the Type B evaluation in Immler et al. (2010, described in Sect. 4.5), is necessary to know if the results obtained in the first method can be regarded as being significantly different.

4.1 Screening with quality assessment

Prior to statistical comparison, irregular or inappropriate data for comparison should be excluded. RS92-GDP, processed by the GRUAN LC, has been checked with a quality assessment algorithm by the GRUAN LC. As the algorithm for iMS-100-GDP in this study is designed with reference to that for RS92-GDP, the algorithm for RS92-GDP is described before that for iMS-100-GDP.

For RS92-GDP, the first quality screening is based on differences between the radiosonde and reference sensor in the ground check. The thresholds are 1.5 hPa for pressure, 1.0 K for temperature, and 1.5 % RH for RH (Sommer, 2013). Data with larger difference are not used. The second screening is based on the quantified uncertainty estimates (Sommer, 2013). This screening is based on the idea that data with uncertainties exceeding the criteria are of questionable reliability and need to be verified individually. The thresholds for uncertainty, $u_{\mathrm{spec}}^{\mathrm{RS}\mathrm{92}}$, are calculated as follows:

$$\begin{array}{}\text{(13)}& {\displaystyle }{u}_{\mathrm{limit}}^{\mathrm{RS}\mathrm{92}}\left(T\right)=\mathrm{0.5}\text{(14)}& {\displaystyle }{u}_{\mathrm{limit}}^{\mathrm{RS}\mathrm{92}}\left(U\right)=\mathrm{0.025}U+\mathrm{2.5}\text{(15)}& {\displaystyle }{u}_{\mathrm{limit}}^{\mathrm{RS}\mathrm{92}}\left(P\right)=\mathrm{0.0004}P+\mathrm{0.6}.\end{array}$$

When uncertainties for 95 % or more of data in the entire profile are within u_limit, the profile is approved and otherwise labelled as “Checked”.

For RH, however, the criterion in Eq. (14) is too strict because it is lower than actual typical values, especially for 50 % RH–70 % RH. Thus, the original formula, as per Eq. (16), is used here.

$$\begin{array}{}\text{(16)}& u_{\mathrm{limit}}^{\mathrm{RS}\mathrm{92}}\left(U\right)=-\mathrm{0.000578}{U}^{\mathrm{2}}+\mathrm{0.0925}U+\mathrm{1.457}.\end{array}$$

However, the lower limit of $u_{\mathrm{limit}}^{\mathrm{RS}\mathrm{92}}\left(U\right)$ is set to 2.5 % RH. Figure 4a–c show the profiles of temperature, RH, and u(U) at 00:00 UTC on 9 November 2018. The blue and red line represent the criteria calculated by Eqs. (14) and (16), respectively. Figure 4d and e show whether the data are within (blue) or exceeding (red) the criteria by Eqs. (14) and (16) at each altitude, respectively. Figure 10f and g show the percentage of data within or exceeding in the whole profile and are, thus, equivalent to the rearrangement of Fig. 10d and e. The percentage of data within the criteria, using Eqs. (14) and (16), is 90 % and 98 % in Fig. 10d and e, respectively. Therefore, this case is classified as “Checked” with the LC criteria (i.e. 95 %). However, it is included in the analysis in this study.

For iMS-100-GDP, the first screening is performed at the ground check under room conditions during preparation. The thresholds of the differences from reference sensors for temperature and RH are ±0.5 ^∘C and 7 % RH, respectively. Sensor values exceeding these criteria are not used for observation.

The second screening checks the contamination or changes in the RH sensor specifications associated with icing. In some soundings with iMS-100, abnormal RH profiles, such as the one shown in Fig. 11 (for 12:00 UTC on 22 September 2017), are sometimes observed. For RS92, there is little need to consider the potential freezing of its humidity sensor since it is heated. On the other hand, the very slow RH decreasing after its passage through the supercooled layer (red shading in Fig. 11b) and relatively high RH values in the stratosphere for iMS-100-GDP are not explained by the hysteresis but are probably due to contamination and changes in the RH sensor specifications related to icing or freezing during passage through supercooled droplet clouds. As checking to determine whether the radiosonde has passed through such clouds is impractical from the RH profile, ice-saturated regions (pink shaded layer in Fig. 11b) are considered. However, as not all data associated with ice-saturated clouds' passage are contaminated, the length of the ice-supersaturated region (ISSR) and the RH and water vapour mixing ratio at the top of the upper troposphere and lower stratosphere (UTLS) are used for screening in this study. ISSR determination is based on saturated water vapour for liquid water and ice that is calculated using Hyland and Wexler equation (Hyland and Wexler, 1983). The observed RH is the ratio of water vapour pressure to e_{s_liq} and is limited to $U_{\mathrm{s}\mathrm{\_}\mathrm{ice}}=e_{\mathrm{sat}\mathrm{\_}\mathrm{ice}}/e_{\mathrm{sat}\mathrm{\_}\mathrm{liq}}$ for the ice phase. Accordingly, ISSR is the layer in which U>U_{sat_ice}. In this study, two adjacent ISSRs, with intervals of ≤60 s or where the minimum of U−U_{sat_ice} is more than −10 % RH (i.e. with no dry layer in the interval), are merged and treated as a single ISSR.

Figure 11 Temperature (a) and RH (b) profiles for RS92-GDP (orange) and iMS-100-GDP (blue) at 12:00 UTC on 22 September 2017. The dashed line in panel (a) represents tropopause pressure. The black line in panel (b) shows the RH profile for ice-saturated air with iMS-100-GDP temperature; air is supersaturated for ice if the observed RH exceeds this value. The pink shaded layer (453.2–232.4 hPa) shows the ISSR for iMS-100-GDP.

For the screening of ice-contaminated profiles, the probability of icing, Pr_ice, is derived using logistic regression analysis (e.g. Cox, 1958) after variable selection from the length of ISSR, RH, and the volume mixing ratio at several levels. The 452 routine observations (twice per day) with a single payload taken from April to November 2018 at Tateno are used as training data. The Python scikit-learn package (Pedregosa et al., 2011) was used for the actual coefficient derivation. The Pr_ice is calculated with the following:

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& {\mathrm{Pr}}_{\mathrm{ice}}=\\ & {\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\exp \left(-\mathrm{22.38}+\mathrm{7.18}\times {\mathrm{10}}^{-\mathrm{3}}{T}_{\mathrm{ISSR}\mathrm{\_}\max }+\mathrm{0.517}{U}_{\mathrm{ST}\mathrm{1}}+\mathrm{0.977}{U}_{\mathrm{ST}\mathrm{2}}+\mathrm{0.105}{W}_{\mathrm{ST}\mathrm{1}}+\mathrm{0.427}{W}_{\mathrm{ST}\mathrm{2}}\right)}},\end{array}\end{array}$$

where subscripts ST1 and ST2 represent the data points at 2000 and 4000 m above the tropopause in geopotential height, U (% RH) is RH, W (parts per million by volume, ppmv) is the water vapour volume mixing ratio, and T_{ISSR_max} is the pass-through period of maximum ISSR in seconds. In this study, profiles with Pr_ice>0.5 are considered to be potentially ice-contaminated data and excluded from the RH analysis below.

The third screening is based on uncertainties for temperature, RH, and pressure as with RS92. The thresholds for these values are derived from Eqs. (18)–(20), respectively. Coefficients are determined empirically.

$$\begin{array}{}\text{(18)}& u_{\mathrm{limit}}^{\mathrm{iMS}}\left(T\right)=\mathrm{0.7}.\end{array}$$

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& u_{\mathrm{limit}}^{\mathrm{iMS}}\left(U\right)=\\ & \left\{\begin{array}{l}-\mathrm{1.438}\times {\mathrm{10}}^{-\mathrm{6}}{U}^{\mathrm{3}}+\mathrm{9.867}\times {\mathrm{10}}^{-\mathrm{4}}{U}^{\mathrm{2}}-\mathrm{1.020}\times {\mathrm{10}}^{-\mathrm{2}}U+\mathrm{3.081}\left(\mathrm{for}\mathrm{daytime}\right)\\ -\mathrm{1.195}\times {\mathrm{10}}^{-\mathrm{5}}{U}^{\mathrm{3}}+\mathrm{2.631}\times {\mathrm{10}}^{-\mathrm{3}}{U}^{\mathrm{2}}-\mathrm{7.308}\times {\mathrm{10}}^{-\mathrm{2}}U+\mathrm{3.351}\left(\mathrm{for}\mathrm{nighttime}\right)\end{array}.\right.\end{array}\end{array}$$

However, the lower limit of $u_{\mathrm{limit}}^{\mathrm{iMS}}\left(U\right)$ is set to 3.3 % RH for both cases.

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{u}_{\mathrm{limit}}^{\mathrm{iMS}}\left(P\right)& =\mathrm{1.609}\times {\mathrm{10}}^{-\mathrm{9}}{P}^{\mathrm{3}}-\mathrm{4.589}\times {\mathrm{10}}^{-\mathrm{6}}{P}^{\mathrm{2}}\\ & +\mathrm{4.120}\times {\mathrm{10}}^{-\mathrm{3}}P+\mathrm{3.810}\times {\mathrm{10}}^{-\mathrm{2}}.\end{array}\end{array}$$

For example, the criteria for temperature and RH uncertainty are shown as black lines in Figs. 4b and d, respectively. In this case, the RH uncertainty exceeds the criteria near 370 and 250 hPa. The ratio criterion is set to 90 %; thus, if more than 90 % of the whole profile has data whose uncertainty does not exceed the threshold, the iMS-100-GDP is included in the following analysis. For the case in Fig. 4, the ratio of data for which the uncertainty is within the threshold is greater than 90 %. On the other hand, Fig. 12 (at 00:00 UTC on 6 December 2019) shows the rejected case with only 84 % of the data within the RH criteria due to the large uncertainty above the 15 hPa level. For the case in Fig. 13 (at 00:00 UTC on 24 November 2017), the large pressure uncertainty layer exceeding the criteria is found between 540 and 230 hPa. Only 82 % of the data are within the pressure criteria, and thus this case is rejected. In this case, the large correlated uncertainty in the middle troposphere is caused by a large vertical component of the dilution of precision (VDOP), i.e. the precision of GPS positioning (not shown). It should be noted that this quality assessment and screening are for this study only and are not authorised as the standard for GDPs. The standard method for quality assessment of GDPs is currently under discussion by the quality task force in the GRUAN community.

Figure 12 As per Fig. 4 but for the case at 00:00 UTC on 6 December 2019.

Figure 13 As per Fig. 4 but for the case at 00:00 UTC on 24 November 2017.

As a result of the screening process described above, the 59 dual sounding flights (29 daytime and 30 nighttime) shown in Table 6) are used as a suitable data set for comparison. Seasonal profiles for temperature, RH, and wind speed are shown in Figs. 15–17, with the seasonal classification (March–April, MA, May–July, MJJ, August–November, ASON, and December–February, DJF) described in Sect. 4.4. The major factors associated with data screening are illustrated in Fig. 14 for each season and daytime/nighttime. The rejection rates are low in DJF but high in MA and MJJ. The primary reasons are the large uncertainties in humidity or pressure, the icing of the humidity sensor for iMS-100, and the bad RS92 sensors, i.e. humidity differences exceeding thresholds to the reference sensor in the ground check. Thus, the data are rejected mainly due to humidity sensor issues. It maybe is related to the climatology at Tateno, where it tends to be dry in winter and humid and highly variable with altitude in MA and MJJ (Fig. 16). The highly variable humidity may cause large uncertainties for the time lag correction or hysteresis correction. Passing through ISSRs increases the likelihood of the humidity sensor icing. For RS92, the ground check of the humidity sensor is performed with the assumption that the humidity of the chamber filled with molecular sieves is 0 % RH, but the molecular sieves in the chamber gradually becomes damp because they take in water vapour from the atmosphere. At Tateno, molecular sieves were replaced when the humidity difference between RS92 sensors, which exceeded 2 % RH and which may have exceeded the RS92-GDP screening criteria.

Table 6 List of dual sounding events with iMS-100 and RS92 used for comparison. The Weather code is according to the World Meteorological Organization (WMO) code table 4677.

^* Truncated due to excessive temperature gaps. RS92 achieved 34 434.6 m (6.1 hPa).

Figure 14 Major factors in the data screening. The term “Adopted” is used for comparison.

Figure 15 Seasonal profiles of iMS-100-GDP. Grey lines represent profiles for individual flights, red represents profile means, and black represents profile medians. Shading represents 25th–75th percentiles.

Figure 16 As per Fig. 15 but for RH.

Figure 17 As per Fig. 15 but for wind speed.

4.2 Timestamp adjustment

Radiosonde observations have timestamps from the relevant sounding system (for iMS-100, this is based on received GPS clock data). As there may be minor discrepancies in balloon launch timestamps, these data are time-adjusted using the temperature profile. In this study, the shift registration in functional data (Ramsay and Silverman, 2005) is used for the adjustment using the temperature data from between 5 and 20 min of the launch time of RS92. Temperature data for each radiosonde, T_i(t), in this period are converted to functional data, x_i(t). Shift values, δ_i, are calculated to minimise the registered curve sum of square error (REGSEE), which is defined as follows:

$$\begin{array}{}\text{(21)}& \mathrm{REGSEE}=\sum _i\int {\left[x_i\left(t+{\mathit{\delta }}_i\right)-\widehat{\mathit{\mu }}\left(t\right)\right]}^{\mathrm{2}}\mathrm{d}t,\end{array}$$

where $\widehat{\mathit{\mu }}\left(t\right)$ is the mean function of x_i(t). For the calculation of δ_i, the scikit-fda package (Carreño et al., 2020) for Python is used, and actual adjustment values (seconds to shift) are derived with the following:

$$\begin{array}{}\text{(22)}& \mathit{\delta }t=\mathit{\delta }{t}_{\mathrm{RS}\mathrm{92}}-\mathit{\delta }{t}_{\mathrm{iMS}}.\end{array}$$

When δt>0, the timestamps for iMS-100 are shifted backward by $\left|\mathit{\delta }t\right|$, and vice versa.

4.3 Data pretreatment

Profiles with >10 % of abnormal data points in the whole profile are excluded via screening, as described in Sect. 4.1, while abnormal data are also seen at individual points in the overall normal profiles (e.g.  superadiabatic lapse rates). Such data points should be excluded from statistical comparison. In this study, superadiabatic lapse rate layers and abnormal wind data immediately before balloon burst are pre-treated for masking.

4.4 Statistical comparison for binned layers based on pressure

After timestamp adjustment, the per second differences between iMS-100 and RS92 measurements were calculated, and the differences were allocated to the 13 pressure layers based on RS92 pressure data ($P_i^{\mathrm{RS}\mathrm{92}}$, where i indicates the time step) based on Kobayashi et al. (2012, 2019). The bins for the 13 layers are listed in Table 7.

Table 7 Pressure range for allocation of iMS-100–RS92 differences based on RS92 pressure (i.e. $\mathrm{bottom}\ge P_i^{\mathrm{RS}\mathrm{92}}>\mathrm{top}$).

Table 8 Terminology for pair checking in independent measurements with identical quantities for consistency.

$A_{i,j}^{\mathrm{iMS}}$ and $A_{i,j}^{\mathrm{RS}\mathrm{92}}$ are values at the time step i for iMS-100 and RS92 in the jth dual sounding data set ($j=\mathrm{1},\mathrm{\dots },M$; M is the number of data sets, here 59), respectively. The differences between the two radiosonde types, $\mathrm{\Delta }{A}_{i,j}=A_{i,j}^{\mathrm{iMS}}-A_{i,j}^{\mathrm{RS}\mathrm{92}}$, are averaged for each pressure layer with the following:

$$\begin{array}{}\text{(23)}& \stackrel{\mathrm{‾}}{\mathrm{\Delta }{A}_{k,j}}={\displaystyle \frac{{\sum }_{i=i_{\mathrm{b},j}}^{i_{\mathrm{t},j}}\mathrm{\Delta }{A}_{i,j}}{N_{k,j}}},\end{array}$$

where k is the layer number ($k=\mathrm{1},\mathrm{\dots },\mathrm{13}$), N_k,j is the number of data points in the layer k, and i_b,j and i_t,j are the bottom and top time step at the jth data set, respectively. The ensemble mean of the difference in each pressure layer is calculated as follows:

$$\begin{array}{}\text{(24)}& \langle \mathrm{\Delta }{A}_k\rangle ={\displaystyle \frac{{\sum }_{j=\mathrm{1}}^M\stackrel{\mathrm{‾}}{\mathrm{\Delta }{A}_{k,j}}}{M}}.\end{array}$$

The ensemble standard deviation of the mean difference for individual pressure layers is as follows:

$$\begin{array}{}\text{(25)}& \mathit{\sigma }\left(A_k\right)=\sqrt{{\displaystyle \frac{{\sum }_{j=\mathrm{1}}^M(\stackrel{\mathrm{‾}}{\mathrm{\Delta }{A}_{k,j}}-\langle \mathrm{\Delta }{A}_k\rangle {)}^{\mathrm{2}}}{M}}}.\end{array}$$

A comparison of the wind data is performed with the indices of the wind speed bias, BIAS, and the root mean square error of the vector difference, RMSVD (CGMS, 2003), for each pressure layer at each dual sounding. BIAS, vector difference (VD), mean vector difference (MVD), standard deviation vector difference (SD), and RMSVD are calculated as follows:

$$\begin{array}{}\text{(26)}& {\mathrm{BIAS}}_{k,j}=\sum {\displaystyle \frac{\mathrm{1}}{N}}(V_{i,j}^{\mathrm{iMS}}-V_{i,j}^{\mathrm{RS}\mathrm{92}})\end{array}$$

$$\begin{array}{}\text{(27)}& {\mathrm{VD}}_{i,j}=\sqrt{(u_{i,j}^{\mathrm{iMS}}-u_{i,j}^{\mathrm{RS}\mathrm{92}}{)}^{\mathrm{2}}+(v_{i,j}^{\mathrm{iMS}}-v_{i,j}^{\mathrm{RS}\mathrm{92}}{)}^{\mathrm{2}}}\end{array}$$

$$\begin{array}{}\text{(28)}& {\mathrm{MVD}}_{k,j}={\displaystyle \frac{\mathrm{1}}{N_{k,j}}}\sum {\mathrm{VD}}_{i,j}\end{array}$$

$$\begin{array}{}\text{(29)}& {\mathrm{SD}}_{k,j}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N_{k,j}}}\sum ({\mathrm{VD}}_{i,j}-{\mathrm{MVD}}_{k,j}}{)}^{\mathrm{2}}\end{array}$$

$$\begin{array}{}\text{(30)}& {\mathrm{RMSVD}}_{k,j}=\sqrt{{\mathrm{MVD}}_{k,j}^{\mathrm{2}}+{\mathrm{SD}}_{k,j}^{\mathrm{2}}},\end{array}$$

where V, u, and v are the wind speed, zonal wind speed, and meridional wind speed, respectively. The mean difference in wind direction (Szantai et al., 2007), $\stackrel{\mathrm{‾}}{\mathrm{\Delta }\mathrm{DIR}}$, is also used.

$$\begin{array}{}\text{(31)}& \stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathrm{DIR}}_{k,j}}={\displaystyle \frac{\mathrm{180}}{N\mathit{\pi }}}\sum \arccos \left({\displaystyle \frac{{\mathit{V}}_{i,j}^{\mathrm{iMS}}\cdot {\mathit{V}}_{i,j}^{\mathrm{RS}\mathrm{92}}}{V_{i,j}^{\mathrm{iMS}}{V}_{i,j}^{\mathrm{RS}\mathrm{92}}}}\right),\end{array}$$

where V is the wind vector. The ensemble mean of these parameters is then calculated as for Eq. (24).

Figure 18 (a) Profiles of ensemble mean temperature differences (iMS-100-GDP minus RS92-GDP) with standard deviations (error bars) for all seasons combined. Red, blue, and black lines represent daytime, nighttime, and all data, respectively. Points are vertically shifted for ease of viewing. Panel (b) is as per panel (a) but represents seasonal profiles for daytime data. The green, red, orange, blue, and black lines represent MA, MJJ, ASON, DJF, and all seasons, respectively. Panel (c) is as per panel (b) but for nighttime data.

Statistics for each pressure layer are calculated for daytime, nighttime, and individual seasons. Due to the safety consideration described in Sect. 3, few dual soundings are performed in July and August, and the weather conditions for August flights are categorised as for autumn rather than summer. Thus, in this study, the flights from August to November are categorised as those for autumn (here, ASON), and flights in previous studies for the seasonal comparison campaigns of radiosondes at Tateno (Kobayashi, 2015; Kobayashi and Hoshino, 2018) are categorised as follows: the summer covered the period from May to June, flights in March and April are categorised with spring (here, MA), those from May to July are categorised with summer (MJJ), and those from December to February are categorised with winter (DJF).

Figure 19 Temperature consistency rank for individual pressure layers.

Figure 20 Distribution of temperature differences for daytime observations. aX is the histogram, bX is the box plot, and cX is the quantile–quantile plot. X indicates the layer for (1) L09 (50–30 hPa), (2) L07 (100–70 hPa), (3) L05 (200–150 hPa), and (4) L03 (500–300 hPa).

4.5 Method for verification of consistency with uncertainties

The uncertainty estimation for RS92 and iMS-100 GDPs are described in Dirksen et al. (2014) and Hoshino et al. (2022a), respectively.

Immler et al. (2010) proposed the terminology for comparing pairs of independent measurements of the same quantity for consistency using estimated uncertainties, as described in the following. Consider two independent measurements, m₁ and m₂, of the same measurand with standard uncertainties, u₁ and u₂, respectively. Under the assumptions that the measurements m₁ and m₂ have independent errors with measurement uncertainties u₁ and u₂, and the difference m₁−m₂ has a normal distribution, the probability that

$$\begin{array}{}\text{(32)}& \left|m_{\mathrm{1}}-m_{\mathrm{2}}\right|>k\sqrt{u_{\mathrm{1}}^{\mathrm{2}}+u_{\mathrm{2}}^{\mathrm{2}}}\end{array}$$

occurs only by chance is roughly 4.5 % for k=2 and 0.27% for k=3. If Eq. (32) is true for k=2, then it is very likely that the two measurements did in fact not measure the same thing, probably due to an unrecognised or unaccounted-for systematic effect in either one or both measurements. Immler et al. (2010) proposed an expression for the degree of consistency, as shown in Table 8. This approach is a Type B evaluation of uncertainty. For the Type A evaluation, Immler et al. (2010) concluded that it is not expected to play an important role within GRUAN, and thus it is not considered in this study.

For a statistical consistency check, the total consistency ranks shown in Table 8 (k<1 is consistent, $\mathrm{1}\le k<\mathrm{2}$ is in agreement, $\mathrm{2}\le k<\mathrm{3}$ is significantly different, or k≥3 is inconsistent) between RS92 and iMS-100 within a specific pressure layer for a particular parameter are estimated as being the 95 % percentile value of consistency ranking numbers within the layer.

5 Results

5.1 Temperature

Figure 18 shows the ensemble mean (lines) and standard deviation (error bar) of temperature differences for daytime and nighttime for all seasons (Fig. 18a), seasonal daytime (Fig. 18b), and seasonal nighttime (Fig. 18c). In the stratosphere, the iMS-100-GDP values are around −0.5 K lower than RS92-GDP for the daytime. For the nighttime, differences are around −0.1 K below the 10 hPa level. Seasonal differences are small.

Figure 19 shows the percentage of the consistency rank in each layer for daytime (Fig. 19a) and nighttime (Fig. 19b). The percentages of significantly different or inconsistent data exceed 50 % in L07 (100–70 hPa) and 30 % in L08 (70–50 hPa) and L09 (50–30 hPa) for daytime observations (Fig. 19a). However, 80 % of the data are in agreement for nighttime comparison (Fig. 19b). Figure 20 shows the distribution of temperature at L09 (50–30 hPa), L07 (100–70 hPa), L05 (200–150 hPa), and L03 (500–300 hPa). Figure 20a1, a2, and a3 show that, in general, differences in L09, L07, and L05 are normally distributed with a sufficient sample size. Accordingly, daytime temperature differences in the stratosphere and the upper troposphere are associated with systematic effects. These differences (especially between 100 and 30 hPa) are attributed to the difference in the solar radiation correction models. Further discussion about the contributions of different radiation heating correction methods to the temperature difference needs other observation-data-like satellites or other types of radiosondes, like RS41. But GNSS-RO-based temperature data are very limited, and no comparative observations have been made at Tateno between the three sondes (iMS-100, RS92, and RS41). Therefore, additional discussion is expected after the results of the comparisons between iMS-100 vs. RS41 and RS92 vs. RS41 are published. However, some samples show significant differences ($>\pm \mathrm{0.5}$K) even in the troposphere or for nighttime soundings (not shown), which are associated either with issues during flights or calibration problems.

Kobayashi et al. (2019) found that the RS-11G-GDP.1 temperature data are about −0.4 K lower than RS92-GDP.2 data for daytime observations in the stratosphere. This means that the iMS-100-GDP temperatures show larger differences from RS92-GDP.2 than the RS-11G-GDP.1 temperatures do. On the other hand, the ratio of data that are evaluated as consistent or in agreement with RS92-GDP.2 temperatures is greater for iMS-100-GDP than RS-11G-GDP.1. This is probably due to the newly included correction for the sensor orientation effects in iMS-100-GDP, which will increase the uncorrelated uncertainty. Further investigation is needed, using intercomparison results between RS-11G and iMS-100 radiosondes.

5.2 Relative humidity

Figure 21 shows the ensemble mean (lines) and standard deviation (error bars) of RH differences for daytime and nighttime for all seasons (Fig. 21a), seasonal daytime (Fig. 21b), and seasonal nighttime (Fig. 21c). This shows that iMS-100 RH is around 1 % RH–2 % RH larger than the RS92 RH around the tropopause and −1 % RH smaller in the stratosphere. Unlike temperature, systematic differences (i.e. biases) between daytime and nighttime soundings are small, but seasonal variations are large in the upper troposphere and the lower stratosphere. Figure 22 shows the ensemble mean of RH differences for six RH ranks. In the lower troposphere (below 500 hPa level), differences are below 2 % RH and exhibit limited correspondence with RH values. In the middle and upper troposphere (500–100 hPa), the difference increases with altitude for 100 % RH–90 % RH. For data with RH ≤10 % RH, the difference is within ±1 % RH for layers above 500 hPa level; iMS-100 RH is wetter in the troposphere and drier in the stratosphere. For data with RH >90 % RH, the data set is limited above 300 hPa and drier than RS92 figures for the middle troposphere (500–300 hPa).

Figure 21 As per Fig. 18 but for RH.

Figure 22 Mean RH difference for six RH ranks.

Kobayashi et al. (2019) found that RS-11G-GDP.1 RH data show about 2 % RH dry tendencies for conditions with RH >90 % RH and about 1 % RH wet tendencies for conditions with RH ≤10 % RH. These different behaviours between RS-11G-GDP.1 and iMS-100-GDP with respect to RS92-GDP.2 in very dry and very wet conditions make the ΔU profiles different in the lower troposphere (below 700 hPa) and in the stratosphere (above 50 hPa; see Fig. 21 in this study and Fig. 11 in Kobayashi et al., 2019). In particular, the 1 % RH wet bias of RS-11G-GDP.1 with respect to RS92-GDP.2 and the dry bias of iMS-100-GDP in the stratosphere has notable differences. The major reason for these differences could be the inclusion of the hysteresis correction in the iMS-100-GDP. The ratio of data that are evaluated as consistent or in agreement with RS92-GDP.2 RH data in the troposphere are greater for RS-11G-GDP.1. This is also due to the inclusion of the hysteresis correction in the iMS-100-GDP because the RH profiles in the troposphere often show rapid changes, as shown in Fig. 16. Further investigation on the differences between RS-11G and iMS-100, resulting from intercomparison flights, is a future task.

Figure 23 shows the percentage of consistency rank in each layer for daytime (Fig. 23a) and nighttime (Fig. 23b). In the troposphere and lower stratosphere (L1–L7, below 70 hPa level), around 10 %–20 % of data are significantly different or inconsistent. The RH profiles for individual flights (Fig. 16) show that RH often shows rapid changes. For the flight at 12:00 UTC (21:00 LT) on 19 April 2019 (Fig. 24), iMS-100 shows a slow tendency in relation to the rapid decrease in RH (e.g. at about 330 hPa) compared to RS92. This difference in the response to rapid changes, especially in case from high to low humidity, is considered to be a reason for the inconsistency of 1 s RH values between the two radiosondes. As described in Sect. 2.2.2, the iMS-100 RH sensor has hysteresis with the large time constant, but the RS92 RH sensor is heated, and its hysteresis is negligible. This difference in the characteristics of the RH sensor could cause a large difference, especially in the rapidly decreasing RH case. In the stratosphere (L8–L13, above the 70 hPa level), RH data from iMS-100 and RS92 seem to be almost in agreement.

Figure 23 As per Fig. 19 but for RH.

Figure 24 Comparison between 500 and 150 hPa levels (L3 to L5) at 12:00 UTC (21:00 LT) on 19 April 2019. (a) Temperature. (b) RH. (c) RH difference.

5.3 Pressure

Figure 25 shows the ensemble mean (lines) and standard deviation (error bars) of pressure differences for daytime and nighttime for all season (Fig. 25a), seasonal daytime (Fig. 25b), and seasonal nighttime (Fig. 25c). The absolute value of the ensemble mean difference is less than 0.4 hPa, but there are cases with large differences in the lower troposphere (below 700 hPa level). This may be attributable to the effect of the pressure differences between the RS92 pressure sensor and the barometer used for surface observation. The pressure of iMS-100, with no pressure sensor, is derived from a recursive calculation via the hydrostatic equation so that the surface pressure is equal to that observed using ground-based barometer. Meanwhile, RS92 involves independent pressure sensor usage, meaning that near-surface pressure may differ between GDPs. A histogram of the RS92 GDP surface pressure error for ground-based barometer content (Fig. 26) shows that the difference is not normally distributed around zero. RS92 pressure tends to be slightly higher than the ground-based barometer content, with a difference median of 0.33 hPa, which is greater than the barometer uncertainty (0.06 hPa for k=1). The consistency check for surface pressure between the barometer and RS92 resulted in significantly different or inconsistent results in 20 cases. The effect of this difference decreases with height but is more noticeable near the ground.

Figure 25 As per Fig. 18 but for pressure.

Figure 26 Histogram for the surface pressure difference (RS92 minus barometer values).

For RS-11G-GDP.1, the pressure is about −0.5 hPa lower for daytime observation in the middle of troposphere, but the difference between iMS-100-GDP and RS92-GDP.2 is statistically small in those layers. At this time, the reason for this difference in pressure is not clear.

5.4 Geopotential height

Figure 27 shows the ensemble mean (lines) and standard deviation (error bars) of geopotential height differences for daytime and nighttime for all seasons (Fig. 27a), seasonal daytime (Fig. 27b), and seasonal nighttime (Fig. 27c). The difference in geopotential height is around 2–3 m in the lower and middle troposphere but becomes larger with altitude above 100 hPa level and becomes about 10 m at 20 hPa. This tendency is attributed to the difference in geoid height, as referenced by iMS-100-GDP and RS92-GDP. As described in Sect. 2.2.3, the grid size of the original geoid model used for the iMS-100 GPS module is $\mathrm{10}{}^{\circ }\times \mathrm{10}{}^{\circ }$, which is replaced with a ${\mathrm{5}}^{\prime }\times {\mathrm{5}}^{\prime }$ model for geometric height calculation in GDP processing. The grid size difference causes a significant discrepancy in geoid and geometric height, especially for the northwestern Pacific basin around Japan. Figure 28a and b show the geoid height for iMS-100-GDP and original iMS-100, respectively, Fig. 28c shows differences in geoid height and typical radiosonde track for each season (green for MA, red for MJJ, orange for ASON and blue for DJF), and Fig. 28d shows geoid correction values for typical seasonal tracks. The difference increases with height for all seasons. The grid size of the geoid model for RS92-GDP is unknown, but as geometric height values are used without modification, the geoid model differences may have caused geopotential height differences.

Figure 27 As per Fig. 18 but for geopotential height.

Figure 28 (a) Geoid height for iMS-100-GDP, (b) iMS-100 original values, (c) difference in geoid height and seasonal typical tracks, and (d) correction values along typical tracks.

The geopotential height difference between iMS-100-GDP and RS92-GDP.2 seems to have no noticeable difference with that between RS-11G-GDP.1 and RS92-GDP.2, although the RS-11G-GDP.1 is not a corrected geoid model. This implies that the geoid model used in RS-11G has enough resolution for GDP.

5.5 Wind

Figure 29 compares wind for the wind speed bias (Fig. 29a), RMSVD (Fig. 29b), and wind direction difference for all seasons (Fig. 29c). The difference is small enough, with BIAS from −0.04 to +0.14 m s⁻¹, RMSVD is less than 1.04 m s⁻¹, except for L13 (above 10 hPa), and $\stackrel{\mathrm{‾}}{\mathrm{\Delta }\mathrm{DIR}}$ is less than 6.4^∘.

Figure 29 Seasonal wind difference parameters with colouring, as per Fig. 18b. (a) BIAS. (b) RMSVD. (c) $\stackrel{\mathrm{‾}}{\mathrm{\Delta }\mathrm{DIR}}$.

Similar to the verification for temperature and relative humidity, the consistency checks in each pressure-binned layer for wind speed and direction are conducted (see Fig. 30). The figure shows that more than 20 % of the data for wind speed are significantly different or inconsistent in the lowest layer (below 700 hPa level) and between 100 and 20 hPa levels. For wind direction, more than 20 % of the data in most layers and more than 50 % of the data in the lowest layer and between 100 and 20 hPa levels are significantly different or inconsistent. This result appears to be inconsistent with the verification result, using the mean differences for each layer mentioned above.

Figure 30 Wind consistency rank for individual pressure layers. (a) Wind speed. (b) Wind direction.

To understand this difference, the wind speed difference and the uncertainty profile at 12:00 UTC on 12 July 2019 are shown in Fig. 31 as an example. The inset figures are an enlarged view of L07 (100–70 hPa), where the wind speed was estimated as being significantly different. In Fig. 31a, which shows the wind speed profile, the overall change trends are consistent for both GDPs, but the RS92-GDP (the orange line) looks smoother than the iMS-100-GDP (the blue line). This difference in the smoothing causes the variability in wind speed differences (black line in Fig. 31b), resulting in some of the data exceeding the range of the synthetic uncertainty threshold (the green line in Fig. 31b for k=2). As a result, the ratio of the data that are estimated as consistent or in agreement is less than 95 % in the layer, and the consistency check gives significantly different or inconsistent results. Although the two GDPs are not in agreement in the layer, the differences due to this smoothing are cancelled out when averaged and, therefore, do not appear as statistical errors. Similar significant difference problems due to smoothing are found in the wind direction profiles. Figure 32 shows the wind direction difference and the uncertainty profile for the same sounding, and the inset figures show the enlarged view of L02 (700–500 hPa). Figure 32 suggests that the significant difference in wind direction, caused by the difference in the degree of the smoothing method, is as with the wind speed. Additionally, for the lowest level (below 700 hPa level, especially in the boundary layer), the differences in the positioning performance of the GPS modules immediately after launch, or the possible unstable posture of the payload due to the double pendulum motion, rig rotation, or swaying, may have caused a slight difference in the wind calculation.

Figure 31 (a) Profiles of wind speed with uncertainty (shades) at 12:00 UTC on 12 July 2019. Blue and orange lines represent iMS-100-GDP and RS92-GDP, respectively. (b) Profiles of wind speed difference (iMS-100-GDP minus RS92-GDP; black line) and synthetic uncertainty (blue and green lines for k=1 and k=2, respectively) are shown.

Figure 32 As per Fig. 31 but for wind direction. The differences in wind direction are calculated as the angle between the wind vectors for both GDPs, which is not considered to be clockwise/anticlockwise.

6 Comparison with a frost point hygrometer

Due to the technical limitations of the RH sensors mounted on operational radiosondes in low temperatures and dry conditions, GRUAN requires comparisons of RH data with values from reference instruments. At Tateno, the comparison flights with the radiosondes and cryogenic frost point hygrometer (CFH; Vömel et al., 2007, 2016) have been conducted twice a year since 2015. However, since the R-23 (HFC-23) liquid cryogen material used to cool the mirror of CFH was regulated under the Montreal Protocol, a new frost point hygrometer, Meisei SKYDEW (Sugidachi, 2019), was adopted in 2020 for cooling the mirror with a Peltier element instead. Figure 33 shows results from the comparison conducted at 06:00 UTC (15:00 LT) on 21 October 2020. Although the differences in RH are significant (with ΔU exceeding the extended uncertainty of k=2) around the tropopause (92.2 hPa, 16 983.0 m, and −73.4 ^∘C) and in the lower part of the troposphere (800–400 hPa), more than 80 % of iMS-100-GDP RH data, except for L02 (700–500 hPa) and L07 (100–70 hPa), are consistent or in agreement. In particular, almost all data above 70 hPa are consistent with SKYDEW. In this study, the uncertainty in SKYDEW RH is not implemented in the consistency check. If the uncertainty in SKYDEW RH is estimated, then the result is expected to be more consistent between the RH of iMS-100-GDP and SKYDEW. Figure 33c shows that the RH of iMS-100-GDP and SKYDEW is generally in agreement for the entire profile and consistent for most of the troposphere.

Figure 33 Comparison of RH profiles between iMS-100-GDP and SKYDEW at 06:00 UTC (15:00 LT) on 21 October 2020. (a) Air temperature from iMS-100. The dashed line shows the tropopause height (92.2 hPa; 16 983.0 m in geopotential height). (b) RH from both. The blue and black lines show iMS-100-GDP and SKYDEW, respectively. SKYDEW RH is calculated from SKYDEW frost/dew point temperature and iMS-100 air temperature. (c) RH difference in iMS-100 from SKYDEW (black) and RH uncertainty for iMS-100-GDP (blue/green for k=1 and 2, respectively).

7 Summary

To characterise GDPs for iMS-100 and RS92, data from dual soundings conducted at Tateno from September 2017 to January 2020 are analysed in this study. The iMS-100 temperature is around 0.5 K lower than RS92-GDP for daytime observations in the stratosphere, and over 50 % of data from between 100 and 70 hPa and over 30 % from between 70 and 30 hPa show significant differences from RS92-GDP. For nighttime observations, the difference is around −0.1 K, with over 80 % of data showing agreement both in the troposphere and the stratosphere. The difference for daytime measurements in the stratosphere is attributed to the correction procedures for solar radiation heating and differences in sensor characteristics.

The iMS-100 RH is around 1 % RH–2 % RH higher in the troposphere and 1 % RH lower in the stratosphere than RS92, but both GDPs are generally in agreement in the troposphere and stratosphere. The difference may be larger in places where rapid RH change occurs. A comparison flight with the SKYDEW and frost point hygrometer, shows that iMS-100 RH agrees well with SKYDEW both in the troposphere and stratosphere.

While there are some cases where significant differences for pressure are observed in the lower troposphere (≥700 hPa) in the consistency check, the mean pressure difference is less than 0.4 hPa. The difference in geopotential height is around 2–3 m in the lower and middle troposphere but increases with altitude above 100 hPa level, from 10 m at 20 hPa. This relationship between height and related differences may stem from differences in the geoid model used for the two GDPs.

In wind comparison, although the consistency check based on uncertainty may estimate the wind speed and direction for both GDPs as being significantly different for each simultaneous data, the comparison parameters show good correspondence; BIAS is between −0.04 and +0.14 m s⁻¹, RMSVD is lower than 1.04 m s⁻¹, except for above 10 hPa, and $\stackrel{\mathrm{‾}}{\mathrm{\Delta }\mathrm{DIR}}$ is smaller than 6.4^∘ in the statistical comparison for each pressure-binned layer. This seemingly contradictory result is mainly due to the difference in the degree of smoothing for the two GDPs.

The modified data processing described in Sects. 2.2.1 and 2.2.2 will be implemented to processing for RS-11G to create new version of RS-11G-GDP, which will be evaluated with intercomparison campaigns written in Kobayashi et al. (2019, RS-11G and RS92) and Kobayashi and Hoshino (2018, iMS-100 and RS-11G). Further direct comparisons between RS-11G-GDP and iMS-100-GDP will be discussed in a future study.

This study involved the evaluation of the characteristics of iMS-100-GDP values with RS92-GDP as a reference, as the latter is certified as a GRUAN data product. GRUAN certification for iMS-100 is underway, and an ongoing analysis of GDP data is considered important for the provision of high-quality products to the climate research/monitoring community. Since February 2020, regular dual soundings of iMS-100 and Vaisala RS41 (the successor to RS92) have been ongoing. RS41-GDP is under development (von Rohden et al., 2022), and the comparison results will be published when sufficient iMS-100-GDP and RS41-GDP data are available. In addition, iMS-100-GDP, RS41-GDP, and M10-GDP (Dupont et al., 2020) are also candidates of the reference data in the WMO 2022 Upper-Air Instrument Intercomparison Campaign (CIMO Task Team on Upper-air Intercomparison, 2020) which will be conducted in 2022. These data will support further evaluation and improvements in iMS-100-GDP.

The interpolation and the estimation of uncertainty for data missing periods are discussed in some articles. For example, Fassò et al. (2020) proposed a method for temperature data, using the Gaussian process, and Colombo and Fassò (2022) attempted to apply it to RH data. These studies will be considered for future improved versions of the iMS-100-GDP.

Appendix A: Abbreviations

EPS	Expanded polystyrene Styrofoam
GBAS	Ground-based augmentation system
GCOS	Global Climate Observing System
GDP	GRUAN data product
GMDB	GRUAN Meta-data Data Base
GNSS	Global Navigation Satellite System
GPS	Global positioning system
GRUAN	GCOS Reference Upper-Air Network
ISSR	Ice-supersaturated region
JMA	Japan Meteorological Agency
JFMT	JMA transmission format
LC	Lead Centre
Meisei	Meisei Electric Co., Ltd.
NCEI	NOAA National Centers for
	Environmental Information
RH	Relative humidity
RMSVD	Root mean square vector difference
SBAS	Satellite-Based Augmentation System
WMO	World Meteorological Organization