High-resolution humidity profiles retrieved from wind profiler radar measurements

The retrieval of humidity profiles from wind profiler radars has already been documented in the past 30 years and is known to be neither as straightforward and nor as robust as the retrieval of the wind velocity. The main constraint to retrieve the humidity profile is the necessity to combine measurements from the wind profiler and additional measurements (such as observations from radiosoundings at a coarser time resolution). Furthermore, the method relies on some assumptions and simplifications that restrict the scope of its application. The first objective of this paper is to identify the obstacles and limitations and solve them, or at least define the field of applicability. To improve the method, we propose using the radar capacity to detect transition levels, such as the top level of the boundary layer, marked by a maximum in the radar reflectivity. This forces the humidity profile from the free troposphere and from the boundary layer to coincide at this level, after an optimization of the calibration coefficients, and reduces the error. The resulting mean bias affecting the specific humidity profile never exceeds 0.25 g kg⁻¹. The second objective is to explore the capability of the algorithm to retrieve the humidity vertical profiles for an operational purpose by comparing the results with observations from a Raman lidar.

1 Introduction

Over the last 30 years, several authors (Gossard et al., 1982, 1998; Tsuda et al., 2001; Bianco et al., 2005; Stankov et al., 1996, 2003; Furumoto et al., 2003, 2007; Klaus et al., 2006; Imura et al., 2007) discussed the possibility of determining the magnitude of the humidity gradient profiles from measurements of zeroth, first and second moments of wind profiler radars (WPR) Doppler spectra in either the ultra high frequency (UHF) or the very high frequency (VHF) range. The method exploits the clear-air scattering properties of electromagnetic waves, depending on the refractive index and consequently on the thermodynamic properties (pressure, temperature, humidity) of clear air. In addition, most of these authors demonstrated the possibility to retrieve humidity profiles by combining radar measurements with simultaneous measurements from other sensors, i.e., in situ radiosonde observations at a poorer time resolution (Tsuda et al., 2001; Stankov et al., 2003; Furumoto et al., 2006) or remote measurements, such as temperature profiles from a radio acoustic sounding system (RASS; Tsuda et al., 2001), observations from a microwave radiometer profiler (Stankov et al., 1996; Gossard et al., 1998; Bianco et al., 2005; Klaus et al., 2006), precipitable water vapor measurements from a global positioning system (GPS) receiver (Gossard et al., 1999; Imura et al., 2007) or a combination of these measurements. The method was considered to be promising for an operational implementation, providing the benefit of a finer time resolution, compared to the conventional radiosonde observations.

However, as far as we know, no successful attempt to apply this method to operational observations has been reported in literature. The first hurdle is related to the lack of self-consistency of the method and the required synergy between various instruments. Another difficulty is represented by the necessity to carry out accurate measurements of the first three radar Doppler spectra moments, with an accurate calibration of the radar-backscattered power (zeroth-order moment) and a careful post-processing of the radar observations. The latter is needed to guarantee that the velocity and width of the spectral peak (first- and second-order moment, respectively) are not disturbed by external contamination such as ground clutter, radio frequency interferences and spurious echoes like birds. In addition, the relationship between the refractive index and the radar reflectivity is valid only in clear-air conditions and this may raise difficulties when rainfall is mixed with clear air¹. In contrast, smaller and saturated particles are considered as cloud particles. Radars emitting around 1.3 GHz (UHF) are more sensitive to rainfall than those emitting around 440 (UHF) or 45 MHz (VHF) and require a careful processing to separate multiple-peak signals in the vicinity of rainfall events. The final hurdle is represented by the fact that the procedure to determine the humidity profile from the humidity vertical gradient measured by the WPR is not straightforward and the different steps of the process are typically a source of accumulated errors.

In this study, we use WPR observations at 1274 MHz (UHF) collected under a variety of atmospheric conditions at mid-latitudes locations and discuss how we cope with the above-mentioned difficulties. UHF measurements are combined with consecutive radiosonde observations, spaced 6 h or longer, to retrieve high-resolution atmospheric humidity profiles. Section 2 provides the theoretical background for the algorithm used to retrieve the humidity profile, which is based on the method developed by Tsuda et al. (2001). In our modified approach, we introduce a new constraint to integrate the humidity gradient vertical profile, which is the level of transition identified by a reflectivity maximum. This level corresponds either to the top of the mixed or residual layer or to a temperature and/or humidity discontinuity. This approach modification was revealed to be crucial to improve the quality of the results. Section 3 illustrates the experimental sites, the encountered meteorological conditions, the instruments involved in the study and the data processing applied on the WPR measurements. Section 4 illustrates the results of our algorithm under three categories of atmospheric conditions (encountered in the three different datasets) and discusses the calibration coefficients used for each dataset in relation to the coefficient values found in the literature. Finally, in Sect. 5, we apply our method to the retrieving of high-resolution humidity profiles between two consecutive radiosonde observations and discuss the possibility of using it as an operational product. For one specific dataset we took advantage of the measurements of a ground-based Raman lidar located close to the WPR to validate the high-frequency (15 min) humidity profiles retrieved by our algorithm.

2 Theoretical background

The refractive index of the air, n, or refractivity $N=(n-\mathrm{1}){\mathrm{10}}^{\mathrm{6}}$, depends on atmospheric thermodynamic properties (Gossard and Sengupta, 1988):

$$\begin{array}{}\text{(1)}& N=\mathrm{77.6}{\displaystyle \frac{P}{T}}+\mathrm{3.73}\times {\mathrm{10}}^{\mathrm{5}}{\displaystyle \frac{e}{T^{\mathrm{2}}}},\end{array}$$

where T is temperature (K), P is atmospheric pressure (hPa) and e is water vapor partial pressure (hPa). We can also express N in terms of specific humidity, q (kilograms of water vapor per kilogram of moist air). q is the parameter we aim to retrieve in the present work. Using the approximation $q=\mathrm{0.622}\frac{e}{P-\mathrm{0.378}e}\simeq \mathrm{0.622}\frac{e}{P}$, Eq. (1) becomes

$$\begin{array}{}\text{(2)}& N=\mathrm{77.6}{\displaystyle \frac{P}{T}}+\mathrm{5.99}\times {\mathrm{10}}^{\mathrm{5}}{\displaystyle \frac{Pq}{T^{\mathrm{2}}}}.\end{array}$$

The vertical gradient of refractivity M (m⁻¹) can be calculated by applying the following linearized equation for small perturbations (Gossard and Sengupta, 1988):

$$\begin{array}{}\text{(3)}& M={\displaystyle \frac{\partial N}{\partial T}}{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}z}}+{\displaystyle \frac{\partial N}{\partial q}}{\displaystyle \frac{\mathrm{d}q}{\mathrm{d}z}},\end{array}$$

where

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle \frac{\partial N}{\partial T}}=-\mathrm{77.6}{\displaystyle \frac{P}{T^{\mathrm{2}}}}-\mathrm{1.2}\times {\mathrm{10}}^{\mathrm{6}}{\displaystyle \frac{Pq}{T^{\mathrm{3}}}},\text{(5)}& {\displaystyle }& {\displaystyle \frac{\partial N}{\partial q}}=\mathrm{5.99}\times {\mathrm{10}}^{\mathrm{5}}{\displaystyle \frac{P}{T^{\mathrm{2}}}}.\end{array}$$

Equation (3) leads to

$$\begin{array}{ll}{\displaystyle }M={\displaystyle \frac{\mathrm{d}N}{\mathrm{d}z}}& {\displaystyle }=-\mathrm{77.6}{\displaystyle \frac{P}{T^{\mathrm{2}}}}{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}z}}-\mathrm{1},\mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}{\displaystyle \frac{Pq}{T^{\mathrm{3}}}}{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}z}}\\ \text{(6)}& {\displaystyle }& {\displaystyle }+\mathrm{5.99}\times {\mathrm{10}}^{\mathrm{5}}{\displaystyle \frac{P}{T^{\mathrm{2}}}}{\displaystyle \frac{\mathrm{d}q}{\mathrm{d}z}},\end{array}$$

which shows that the refractivity vertical gradient consists of 3 different terms: the temperature gradient term (term 1), the humidity term (term 2) and the humidity gradient term (term 3). As underlined by Tsuda et al. (2001), term 1 may be dominant under dry conditions such as winter conditions at mid-latitudes or in the upper atmosphere. They also found that term 2 usually contributes less than 10 %, even close to the surface where q is the larger. The dominant contribution of term 3, depending on $\frac{\mathrm{d}q}{\mathrm{d}z}$, allows one to solve very easily the first-order differential Eq. (4) (Tsuda et al., 2001). Some authors (Gossard et al., 1998; Stankov et al., 2003) assume that the partial derivatives in Eqs. (4) and (4) are constant and can be estimated from standard atmosphere profiles. This assumption also imposes the neglect of the contribution by the second term. In this paper we do not make this assumption, especially because one of our datasets is characterized by relatively moist conditions near the surface, which forces us to consider all three terms in Eq. (4). As also demonstrated by Tsuda et al. (2001), the differential Eq. (4) can be solved after introducing the Brunt–Väisälä frequency N_BV (s⁻¹) given by

$$\begin{array}{}\text{(7)}& {\left(N_{\mathrm{BV}}\right)}^{\mathrm{2}}={\displaystyle \frac{g}{\mathit{\theta }}}{\displaystyle \frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}z}}=g{\displaystyle \frac{\mathrm{d}\ln \mathit{\theta }}{\mathrm{d}z}},\end{array}$$

where $\mathit{\theta }={\left(\frac{\mathrm{1000}}{P}\right)}^{\mathrm{2}/\mathrm{7}}$ and g are the air potential temperature (K) and the acceleration of gravity (9.8 m s⁻²), respectively.

This leads to

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }M=\mathrm{5.99}\times {\mathrm{10}}^{\mathrm{5}}{\displaystyle \frac{P}{T^{\mathrm{2}}}}{\displaystyle \frac{\mathrm{d}q}{\mathrm{d}z}}-\mathrm{1.2}\times {\mathrm{10}}^{\mathrm{6}}{\displaystyle \frac{Pq}{T^{\mathrm{2}}}}{\displaystyle \frac{(N_{\mathrm{BV}}{)}^{\mathrm{2}}}{g}}\\ \text{(8)}& {\displaystyle }& {\displaystyle }-\mathrm{77.6}{\displaystyle \frac{P}{T}}{\displaystyle \frac{(N_{\mathrm{BV}}{)}^{\mathrm{2}}}{g}},\end{array}$$

which can be rewritten (Tsuda et al., 2001) in the form

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\mathrm{d}q}{\mathrm{d}z}}+A\left(z\right)q=B\left(z\right),\end{array}$$

where

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }A\left(z\right)=-\mathrm{2}{\displaystyle \frac{(N_{\mathrm{BV}}{)}^{\mathrm{2}}}{g}},\text{(11)}& {\displaystyle }& {\displaystyle }B\left(z\right)=\mathrm{1.65}{\displaystyle \frac{T^{\mathrm{2}}}{P}}M+{\displaystyle \frac{\mathrm{1}}{\mathrm{7750}}}\left({\displaystyle \frac{\mathrm{d}T}{\mathrm{d}z}}+\mathrm{\Gamma }\right),\end{array}$$

where $\mathrm{\Gamma }=\mathrm{9.8}\times {\mathrm{10}}^{-\mathrm{3}}$ K m⁻¹ is the dry adiabatic temperature lapse rate. Tsuda et al. (2001) provided the following final form for Eq. (5):

$$\begin{array}{}\text{(12)}& q\left(z\right)={\mathit{\theta }}^{\mathrm{2}}\underset{zo}{\overset{z}{\int }}\left(\mathrm{1.67}\times {\mathrm{10}}^{-\mathrm{6}}{\displaystyle \frac{M{T}^{\mathrm{2}}}{P}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{7750}}}{\displaystyle \frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}z}}\right){\mathit{\theta }}^{-\mathrm{2}}\mathrm{d}z+q_o,\end{array}$$

with q_o being the humidity at level z_o, where the integration is initialized.

Following Tsuda et al. (2001), Furumoto et al. (2006), Klaus et al. (2006) and Imura et al. (2007) used Eq. (6) to compute humidity profiles. We are also using this equation in the present work.

The next step is to relate M to the radar characteristics. In clear-air or cloudy conditions (precipitation-free atmosphere), UHF-range and VHF-range profilers detect the fluctuations of refractive index with a scale of one-half the radar wavelength, through the following expression (Ottersten, 1969):

$$\begin{array}{}\text{(13)}& \mathit{\eta }=\mathrm{0.38}{C}_n^{\mathrm{2}}{\mathit{\lambda }}^{-\mathrm{1}/\mathrm{3}},\end{array}$$

where η is the volume reflectivity for the turbulence echo (m⁻¹), depending on the radar return signal power, λ is the wavelength and $C_n^{\mathrm{2}}$ is the turbulence structure parameter (m${}^{-\mathrm{2}/\mathrm{3}}$) for the radar refractive index. Gossard et al. (1982, 1998) found that, for homogeneous isotropic turbulence in a horizontally homogeneous medium with vertical gradients of mean properties, the squared vertical gradient of potential refractivity (potential refractivity is the value of N for an air parcel moving adiabatically from its ambient level to the reference level – 1000 hPa – without loss or gain of moisture) is

$$\begin{array}{}\text{(14)}& {\left({\displaystyle \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}z}}\right)}^{\mathrm{2}}={\left({\displaystyle \frac{L_w}{L_{\mathit{\varphi }}}}\right)}^{\frac{\mathrm{4}}{\mathrm{3}}}{\left({\displaystyle \frac{\mathrm{d}V}{\mathrm{d}z}}\right)}^{\mathrm{2}}{\displaystyle \frac{C_n^{\mathrm{2}}}{C_w^{\mathrm{2}}}}.\end{array}$$

According to Tatarskii (1971), $\frac{L_w}{L_{\mathit{\varphi }}}$ is the ratio of two outer length scales, for shear (L_w) and for potential refractive index (L_ϕ). Gossard et al. (1982) provided an empirical formulation for this ratio which they refined in Gossard et al. (1998). They found this ratio to be very small in stable layers and large in zones with near-neutral stability, with values ranging between 2 and 6. In this expression $\frac{\mathrm{d}V}{\mathrm{d}z}$ is the vertical shear of the horizontal wind vector:

$$\begin{array}{}\text{(15)}& {\left({\displaystyle \frac{\mathrm{d}V}{\mathrm{d}z}}\right)}^{\mathrm{2}}={\left({\displaystyle \frac{\mathrm{d}u}{\mathrm{d}z}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathrm{d}v}{\mathrm{d}z}}\right)}^{\mathrm{2}},\end{array}$$

where u and v are the horizontal components of the wind. In Eq. (8), $C_w^{\mathrm{2}}$ is the structure parameter of vertical velocity, which can be expressed in term of the dissipation rate of the turbulent kinetic energy ϵ (m² s⁻³) through the expression $C_w^{\mathrm{2}}=\frac{\mathrm{4}}{\mathrm{3}}\times \mathrm{2.1}\times {\mathit{ϵ}}^{\frac{\mathrm{2}}{\mathrm{3}}}$ (where 2.1 is the Kolmogorov constant).

Assuming that $M^{\mathrm{2}}={\left(\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}z}\right)}^{\mathrm{2}}$, Eq. (8) gets the form

$$\begin{array}{}\text{(16)}& C_n^{\mathrm{2}}={\mathit{\alpha }}^{\mathrm{2}}{\displaystyle \frac{{\mathit{ϵ}}^{\mathrm{2}/\mathrm{3}}{M}^{\mathrm{2}}}{{\left(\frac{\mathrm{d}V}{\mathrm{d}z}\right)}^{\mathrm{2}}}}.\end{array}$$

Here, the first three moments of the radar spectral data are involved. Specifically, $C_n^{\mathrm{2}}$ is related to the zeroth moment since the volume reflectivity η in Eq. (7) is proportional to the backscattering signal power. The coefficient of proportionality K between η and the backscattering signal power depends on the radar properties, including the antenna efficiency, receiver bandwidth, system noise power and losses in the transmission lines. In most studies, K is not known since the radar is not calibrated, so K is included in the term α² of Eq. (10).

The horizontal wind V is determined from the first-order moment, i.e., the Doppler shift in the spectral data, obtained at the end of the radar computation process. ϵ is related to the estimation of the second-order moment, i.e., the broadening of the radar spectra.

Let us now discuss the coefficient α². Assuming the radar to be calibrated, so that K is known and the structure coefficient $C_n^{\mathrm{2}}$ is an absolute value, Gossard et al. (1998) took a constant value for α² of 0.44 regardless of the M² profile and the considered level in the profile. This value is obtained by comparing refractive index profiles obtained from balloon measurements to those estimated by a 440 MHz UHF radar, after an accurate post-processing of the first- and second-order moments of the spectral data and an accurate radar calibration (that provided K). In fact, they measured an average value for $\frac{L_w}{L_{\mathit{\varphi }}}$ of 4, which, after Gossard et al. (1982), should be considered as a conservative value under steady conditions.

Following Ottersten (1969), Gossard et al. (1982) and Stankov et al. (2003) used another expression for the turbulence structure parameter, based on the energy equation:

$$\begin{array}{}\text{(17)}& C_n^{\mathrm{2}}=a^{\mathrm{2}}{\displaystyle \frac{{\mathit{ϵ}}^{\mathrm{2}/\mathrm{3}}{M}^{\mathrm{2}}}{{\left(\frac{\mathrm{d}V}{\mathrm{d}z}\right)}^{\mathrm{2}}}}(\mathrm{1}-R_{\mathrm{f}}{)}^{-\mathrm{1}},\end{array}$$

where R_f is the flux Richardson number. They indicate that this relation is valid under fairly general assumptions inside regions of large kinetic energy transfer from shear into turbulence. A value of 2.8 is used for a².

For the specific case of a free troposphere which is hydrostatically stable most of the time and where the turbulence is known to be intermittent (in time and space), VanZandt et al. (1978) refined this relation by inserting F, a “filling factor”, which accounts for the turbulent fraction of the backscattering volume. Based on the introduction of this factor, the expression for the turbulence structure parameter has the form

$$\begin{array}{}\text{(18)}& C_n^{\mathrm{2}}=\mathrm{2.8}{\displaystyle \frac{{\mathit{ϵ}}^{\mathrm{2}/\mathrm{3}}{M}^{\mathrm{2}}}{{\left(\frac{\mathrm{d}V}{\mathrm{d}z}\right)}^{\mathrm{2}}}}(\mathrm{1}-R_{\mathrm{f}}{)}^{-\mathrm{1}}{F}^{\mathrm{1}/\mathrm{3}}.\end{array}$$

VanZandt et al. (1978) provided an estimation of F based upon a simple model for the statistical distribution of wind shear and potential temperature gradient. The comparison of Eqs. (10) and (12) reveals that α² should depend on the stability conditions in both the atmospheric boundary layer (ABL) and free troposphere and also on the filling factor in the free troposphere. In agreement with the concept of filling factor, Gossard et al. (1999) noticed that α² could depend on anisotropy or turbulence unsteadiness. Tsuda et al. (2001) considered that α² is not constant and requires calibration, even if the radar calibration coefficient K is known. In the present work, we will pay special attention to the coefficient α².

To summarize, the first three moments of the radar spectrum allow us to determine M² (Eq. 10), provided that α² is known. The sign of M remains, however, ambiguous. In most cases it is negative, but it can become positive under clouds or in locally dry layers. Supposing that the ambiguity is resolved, the computation of q through Eq. (6) requires the knowledge of the initial condition q_o and of the potential temperature and pressure profiles (or at least an estimation of the air density at the surface to derive the pressure from the temperature profile). As underlined by Klaus et al. (2006), the estimation of a single humidity profile requires simultaneous measurements of four data inputs, in addition to the radar measurements.

As already mentioned in the Introduction, several authors proposed to get the temperature profiles from simultaneous RASS observations and the initial density and humidity conditions from a meteorological station at the surface. This, however, raises an issue since the first radar gate is not at the surface and since the signal-to-noise is sometimes disturbed at the lower gates.

Other authors used the combination of radar wind profilers and ground-based microwave radiometers, the latter providing both temperature and humidity profiles, but at a coarser vertical resolution than the profiler vertical resolution. The first advantage of this method is represented by the possibility to calibrate the radar estimate of the $\mid M\mid$ value using the radiometer integrated humidity on the air column. The second advantage is represented by the possibility to determine the sign of M from the radiometer humidity profile, as in fact this sign is the same as the sign of the humidity gradient or negative if the humidity gradient is negligible.

The integrated value of humidity on the air column can also be obtained from the GPS data and is sometimes used as a constraint for the calibration of $\mid M\mid$, provided that at least one in situ humidity profile is available (e.g., from a radiosonde) to calculate the amount of the integrated humidity in the height range where the profiler identifies echoes relative to the total column integrated humidity. This method also assumes that this contribution is constant, which is a rough hypothesis and difficult to fulfill, especially because of the variability in the radar detection height.

Finally, another method consists in using radiosonde profiles with low time resolution (6 or 9 hourly) to initialize q, calibrate M and provide the temperature and pressure profiles. Then the objective with the radar is to provide intermediate profiles with a finer time resolution. This is the method we decided to use in the present work. In addition, we propose two refinements. Usually, the integration shown in Eq. (6) is initialized at the lowest gate measured by the radar, or, if the latter is too high and observation not available, at the upper boundary where q becomes nearly zero (Tsuda et al., 2001). Here, we make two integrations. The first is started near the first radar gate, which is between 150 and 375 m above ground level (a.g.l.) according to the radar and the measurement mode (see next section), while the second is started at the upper boundary of the radar profile, which is never the same since it depends on the backscattered echo power (the detectability is enhanced under moist conditions and reduced under very dry conditions). Both integrated profiles are adjusted (by joining the points) at a characteristic height that we will call Hlim. Hlim may correspond to the top level of the ABL, under unstable conditions, or at least to a strong moisture gradient under neutral or stable conditions. This level is easily detected by the WPRs, since it corresponds to a local maximum in the radar reflectivity profile. Hlim has been used for the past 20 years to monitor the depth of the mixed ABL (Angevine et al., 1994; Heo et al., 2003).

Since this level can demarcate two regions of the low troposphere, potentially characterized by drastically different turbulent conditions, we decided to compute α² independently in both part of the profiles: below and above Hlim. α² was determined in each of these regions by computing $M_{\mathrm{RS}}^{\mathrm{2}}/M_{\mathrm{r}}^{\mathrm{2}}$. $M_{\mathrm{RS}}^{\mathrm{2}}$ is calculated from the radiosonde data through Eq. (4) and $M_{\mathrm{r}}^{\mathrm{2}}$ from the WPR data through Eq. (10). This calibration step could have exempted us from calibrating the WPRs. Finally, it is to be pointed out that, in our method, we also use the sign of the balloon humidity gradient to determine the sign of M.

3 Experiments, sites, instruments and radar data processing

3.1 Field experiments, measurement sites, instruments and radar data processing

Data were collected during three field experiments characterized by drastically different atmospheric conditions. The first experiment, the Boundary-Layer Late Afternoon and Sunset Turbulence (BLLAST; Lothon et al., 2014), took place during the summer period in June 2011, at the Lannemezan Atmospheric Research Center (43.13^∘ N, 0.13^∘ E; elevation 595 m) at the Pyrenean foothills, and was characterized by typical fair weather convective boundary layer conditions found at mid-latitudes. The second campaign was the first Special Observing Period of the Hydrological cycle in Mediterranean Experiment (HyMeX SOP1), which took place over the western part of the Mediterranean basin during fall 2012 (Ducrocq et al., 2014). In this case the radar was installed in an atmospheric “supersite” located in Candillargues (43.6^∘ N, 4.07^∘ E; elevation 2 m) in the proximity of the seaside. The aim of the experiment was to study the upstream dynamical conditions linked to the initiation of strong rainfall events inland. The preferred situations were those when warm, moist and unstable air masses were advected from the sea. The third experiment was the second Special Observing Period of HyMeX (HyMeX SOP2), which was held during winter 2013, to study the mechanisms of air–sea exchanges in case of strong offshore winds (Estournel et al., 2016). The radar site was the same as for HyMeX SOP1.

3.2 Instruments

During these three experiments, intensive radiosounding (RS) operations (3 or 6 hourly) were performed with the purpose to monitor the atmospheric diurnal cycle. MODEM M2K2 and Vaisala RS92 (OYJ DigiCORA V3.64/RS92-SGP) radiosondes were used for BLLAST and HyMeX, respectively. The accuracy of the RS92 humidity measurements has been assessed by Miloshevich et al. (2006), based on the comparison of these measurements with those simultaneously performed by a reference sensor of known absolute accuracy (a cryogenic frost point hygrometer) deployed on the same balloon. In the low troposphere, as considered in the present study, they found that the mean accuracy in the relative humidity measurements (with respect to the absolute sensor measurement) is always lower than 5 %. In 2010, the World Meteorological Organization conducted a new intercomparison experiment in China (Nash et al., 2011). According to these authors, the Vaisala RS92 version tested in China showed a systematic error of less than 2 % and a random error of ≃ 5 % in relative humidity measurements in daytime and nighttime, from the surface to the low stratosphere, in clear-air or cloudy conditions. During this experiment, the MODEM M2K2 radiosonde was also compared to a set of different radiosondes, including the Vaisala RS92. According to Nash et al. (2011), MODEM nighttime measurements had large positive biases (larger than 10 %) for most of the time in the lower and middle troposphere. MODEM radiosondes were also found to suffer from evaporative cooling when emerging from a cloud layer and to overestimate relative humidity by 6 % at night for relative humidities in the range 90–100 %, i.e., in cloudy air. Luckily the BLLAST measurements, in contrast to the HyMeX observations, were systematically carried out in clear air, which mitigated the uncertainties in the MODEM measurements.

WPRs in the framework of BLLAST and HyMeX were deployed by the Centre de Recherches Atmosphériques (Laboratoire d'Aérologie) and the Centre National de Recherches Météorologiques (Météo-France), respectively. Both profilers are five-beam model PCL 1300 manufactured by Degreane. A detailed description of the WPRs, their main working parameters and the processing methods are provided in Saïd et al. (2016). Both radars were operated almost continuously at both sites over long periods (exceeding 18 months). For the purpose of this work, we will concentrate on the data collected from 19 June to 6 July 2011 during BLLAST, from 13 September to 5 November 2012 during HyMeX SOP1 and from 1 February to 15 March 2013 during HyMeX SOP2. Two operational modes were considered. The first mode (low mode), associated with a pulse length of 0.5∕1 µs, sampled the lower troposphere from 75∕150 m a.g.l. to 5∕5.7 km a.g.l. (for BLLAST/HyMeX, respectively). The second mode (high mode) was specified for higher altitude sampling, from 150 m to 8 km a.g.l., and a pulse length of 2.5 µs was considered for both experiments. The vertical resolution was 75∕150 m (for BLLAST/HyMeX, respectively) in low mode and 375 m in high mode. The high mode was oversampled and provided data every 150 m. The inter-pulse period (IPP) was 40∕45 µs (for BLLAST/HyMeX, respectively) in low mode and 80 µs in high mode for both experiments. The radar beam was steered into five directions: one vertical and four oblique directions at zenith angle of azimuths spaced 17 and 90^∘. The beam width was 8.5^∘, which was narrow enough to enable accurate measurements of the Doppler spectral width in the low troposphere.

In the frame of HyMeX SOP1, the ground-based University of Basilicata Raman lidar system BASIL (Di Girolamo et al., 2009) was deployed in Candillargues and operated from 5 September to 5 November 2012, collecting more than 600 h of measurements, distributed over 51 measurement days and 19 intensive observation periods (Di Girolamo et al., 2016). BASIL makes use of a frequency-tripled Nd:YAG laser source, emitting pulses at 355 nm, with a single pulse energy of 500 mJ and a pulse repetition frequency of 20 Hz (average power at 355 nm: 10 W). The receiver consists of a Newtonian telescope (primary mirror diameter: 45 cm, f∕2.1). The major feature of BASIL is represented by its capability to perform high-resolution and accurate measurements of atmospheric temperature and water vapor mixing ratio (kilograms of water vapor per kilogram of dry air), both in daytime and nighttime, based on the application of the rotational and vibrational Raman lidar techniques, respectively (Di Girolamo, 2004; Di Girolamo et al., 2009). Besides temperature and humidity, BASIL can also provide measurements of particle backscatter, extinction and depolarization at several optical wavelengths. Based on an integration time of 5 min and a vertical resolution of 150 m (which are the resolutions of the lidar data used in this paper), the typical daytime precision in water vapor mixing ratio measurements is 0.2 g kg⁻¹ up to 3 km and 0.3 g kg⁻¹ up to 5 km, while the typical nighttime precision is 0.05 g kg⁻¹ up to 3 km and 0.005 g kg⁻¹ at 10 km (Di Girolamo et al., 2016). In the following, we will use the symbol q for both specific humidity and water vapor mixing ratio, since the percentage deviation between the two is rarely exceeding 1 %, even in case of large humidity concentrations, which is far less than the systematic and statistical uncertainties affecting the lidar mixing ratio measurements. During HyMeX SOP1, BASIL measurements of humidity were calibrated based on the comparison with simultaneous radiosondes (the RS mentioned above) launched from a facility located approximately 100 m away from the lidar. A mean calibration coefficient was estimated by comparing BASIL and radiosonde data at all times when BASIL was running (approximately 50 comparisons).

3.3 Radar data processing

The radar time series were coherently averaged to reduce the computing time while preserving signal detectability. The number of coherent integrations (NCI) was calculated for each cycle of 10 beams (five-beam steering in two modes), according to the wind speed measured during the former cycle, to optimize the Nyquist interval. To obtain radial velocity spectra, 128-point Fourier transforms of the finite time series were applied. Finally, an incoherent integration of 30 consecutive spectra was made to improve the signal detectability. In case of rainfall, IPP was enlarged, NCI changed and backscattered power reduced to avoid second trace echoes or saturation. This first step in raw data processing was made with the software provided by the manufacturer and implemented in the radar sites. One full cycle was achieved every 4 to 5 min (according to NCI) and represented 12 to 15 min observations.

The second step, called consensus, was processed in real time and also post-processed at the Laboratoire d'Aérologie. The purpose of the consensus data processing is to determine the meteorological spectral peak among the four more powerful peaks of the Doppler spectra at each range gate. The method is fully described in Saïd et al. (2016). Special care was devoted to the separation between the meteorological peak and ground clutter echoes or the separation of individual echoes from a multiple-echo peak. This was decisive to provide an accurate spectral width of the Doppler peak, used to compute the dissipation rate of turbulent kinetic energy, ϵ. Furthermore, special tests were performed to separate clear-air spectra from precipitation spectra, using information from the four moments of the vertical velocity. Finally, in the data post-processing we manually flagged out the spectra disturbed by birds that were frequent at night during the fall and late winter seasons of HyMeX, which coincided with migration periods.

3.4 Processing of the dissipation rate of turbulent kinetic energy

The dissipation rate of turbulent kinetic energy ϵ was computed according to the method and coefficients proposed by Jacoby-Koaly et al. (2002). Following Doviak and Zrnic (1984), Hocking (1988), Gossard et al. (1998) and White et al. (1999), they determined ϵ through the estimation of the broadening of the Doppler spectrum peak. The broadening of the spectrum also had to be corrected for contributions due to shear, to the antenna beamwidth and to the filtering effects of the Doppler spectrum. We chose to derive ϵ from a combination of the estimations obtained from the vertical velocity spectrum with the median of the estimations obtained from the oblique velocity spectra. This method had been previously assessed by Jacoby-Koaly et al. (2002).

3.5 Radar calibration

Both radar were calibrated according to the method proposed by Campistron and Réchou (2012) and improved by Campistron et al. (2013). The calibration is based on the comparison between rain rate measured by the profiler and rain gauge at the ground. The height of the radar data is taken as low as possible considering signal saturation, receiver linearity and ground clutter. Usually the best level is found around 600 m a.g.l. Long-lasting stratiform precipitation periods were chosen to avoid the presence of strong vertical air velocities. Also, high relative humidity periods were selected to minimize rain modification during its fall.

Following Ulrich (1983), Campistron et al. (2013) assumed that raindrop size distribution follows a gamma function with two parameters that have to be determined. The drop fall speed in still air can be related to the droplets' diameter taking into account the change of density with height (Atlas et al., 1973; Foote and Du Toit, 1969). The parameters of the gamma distribution can be obtained using the mean vertical velocity and the radar reflectivity factor (Chu and Su, 2008). Finally, the rain rate is derived by integrating the droplets' distribution over the diameter interval supposed to extend from zero to infinity.

During each of the rainfall events selected for the calibration, the radar constant was modified until the best agreement was found between rain gauge and radar measurements. An average of the results was retained as the final calibration constant. As stated before, the WPRs used in the framework of BLLAST and HyMeX had provided longer datasets than the datasets we refer to for the present work. That is why we had no difficulty in finding several stratiform conditions to achieve the calibration. The calibration was done for the BLLAST and HyMeX radars in low mode and for the BLLAST radar also in high mode. For each radar and each operation mode, the variability of the calibration coefficients K induced by the calibration method never exceeded 12 %, which is small relative to the variability of the coefficients α² that will be discussed further on (see Sect. 4.3).

3.6 Data conditioning for the humidity gradient retrieval

During BLLAST, the balloon took around 4 min to cross the ABL, whose top level was typically situated at level 2000 m above sea level (a.s.l.; 1400 m a.g.l.). Within the ABL, the wind remained weak (around 4 m s⁻¹), which corresponded to the preferred conditions for the experiment (clear air, anticyclonic conditions). A nocturnal low level jet sometimes occurred at night, but it was seldom stronger than 6 m s⁻¹. During the first 4  min, the balloon drifted horizontally by 1 km from the release site, located nearby the radar site (150 m). The balloon reached the typical maximum height of the BLLAST radar soundings (4 km) in 17 min (from the release), which corresponds approximately to three full cycles of the profiler. Therefore, for the humidity gradient comparison we used 15 min of radar data, during which the balloon drifted horizontally between 5 and 15 km from the radar site, according to the wind conditions in the free troposphere. At a level of 2 km, the horizontal coverage of radar observations was as small as 625 m, while it was twice as large at a height of 4 km.

During HyMeX SOP1, as anticipated before, the lowest atmospheric layers were mainly characterized by marine conditions, with southerlies or easterlies reaching 10 to 20 m s⁻¹ in the lowest 1–2 km (a.g.l. and a.s.l.). At this time of the year (fall season), the conditions were not propitious for sea-breeze development and mixed boundary layers had few opportunities to develop. Above 1 or 2 km, the wind conditions changed to the typical mid-latitude westerlies. The wind speed was never exceeding 25 m s⁻¹ at 5 km, which was the maximum height reached by the radar echoes during SOP1 (moister conditions than during BLLAST and so better detectability of the echoes). The balloons typically took 5.5 and 14.5 min to reach the heights of 2 and 5 km, respectively. 14.5 min correspond to almost three full cycles of the radar. Due to the shear in the wind direction, the maximum drift of the balloon during its ascent was 10 km from its release point.

The radiosonde data were averaged by slices of 75 or 150 m (according to the radar vertical resolution), centered on the radar gates levels. This was revealed to be a better choice than interpolating the radar data but could yield some slight discrepancies in case of sharp humidity gradients.

We carefully interpolated the radar measurements to fill some gaps in the data. Profiles with gaps larger than 750 m were excluded from further analysis. Second-order radar estimations (vertical shear and ϵ) were smoothed to avoid sharp local derivatives. Since 1274 MHz WPRs are sensitive to both turbulence and raindrops, individual profiles were also checked to remove those contaminated by precipitation echoes (essentially during HyMeX SOP1). However, we kept the profiles for which the precipitation echoes were confined within a limited number of levels and the measured vertical velocity did not exceed 0.5 m s⁻¹.

Figure 1 Vertical profiles of refractivity gradient (a, b) and humidity (c, d) using a systematic negative sign for M (a, c) or after assigning to M the sign of M provided by the RS observations (b, d). RS values (q-RS) are red solid lines and radar values (q-Radar) are black solid lines. The thin red lines (int-q-RS) identify the humidity profiles retrieved from the integration of M_RS after averaging the RS observations by slices of 75 m, to match the vertical resolution of the radar. The red dashed line (qsat-RS) is for the saturated humidity profile. The horizontal dashed line delineates Hlim, the transition level used to separate the upper and lower part of the profile. The two values of α² (calibration coefficients), over or below Hlim, are also indicated.

4 Radar humidity profiles versus radiosonde profiles

4.1 Some adjustments to improve the method

Before providing detailed humidity profiles between consecutive radiosoundings, we had to check how the humidity profiles retrieved by the radar were consistent with the initial radiosonde profiles. This took the longest time to perform since it required several adjustments.

The first adjustment was already mentioned in Sect. 2. We solved the sign ambiguity on M (since the radar provides M²) by assigning to M the sign provided by the RS observations (Eq. 4). Figure 1 illustrates M and q profiles before and after the correction. We had initially assigned a negative sign to M, which had revealed to be relevant for most of HyMeX SOP1 conditions, where the source of humidity is close to the surface (maritime air masses) and decreases with the height. There were, however, situations, during either BLLAST or HyMeX, when drier layers disrupted this negative gradient. We show an example from BLLAST in Fig. 1. In Fig. 1a, M is computed from Eq. (4) for the RS observations (red line) and extracted from Eq. (10) for the radar (black line), assuming a constant negative sign on the whole profile. The resulting radar humidity profile shown in Fig. 1c rapidly deviates from the observed profile (RS profile) and maintains the deviation, in both the upper or lower parts of the profile. The downward and upward integrations computed to retrieve q are connected at Hlim = 1745 m, which corresponds to a maximum in the radar reflectivity profile and to a change in the observed humidity profile. The approach to select this level will be discussed later on. The successive changes in the sign of M are taken into account in Fig. 1b. The resulting radar humidity profile is clearly improved (Fig. 1d). The thin red lines in Fig. 1c and d represent the humidity retrieval as obtained from the integration of M_RS, after averaging the RS observations by slices of 75 m, to match the vertical resolution of the radar. The discrepancy between the two red lines can be considered as a systematic error associated with the loss in vertical resolution linked to the 75 m averaging.

Occasional negative values of q are put to zero, especially during BLLAST when very dry layers were observed. Similarly, some unexpected large values of q were put to the humidity saturated value provided by the RS. This occurred during HyMeX SOP1 when moist conditions were frequently encountered. This limitation enabled us to minimize the error accumulation in case of divergence of the integration. We will provide illustrations of such situations in the following.

Figure 2 Radar turbulence structure parameter $C_n^{\mathrm{2}}$ (a, d) and humidity profiles (same details as in Fig. 1c) for different Hlim levels. Panels (b) and (e) illustrate the WPR specific humidity profile (q-Radar) obtained considering the value of Hlim corresponding to the dominant $C_n^{\mathrm{2}}$ peak observed in panels (a) and (d), respectively. Panels (c) and (f) illustrate the WPR specific humidity profile obtained considering the value of Hlim corresponding to a different relative maximum of the $C_n^{\mathrm{2}}$ profile. In panels (b), (c), (e) and (f), the RS specific humidity profile (q-RS), the saturation specific humidity profile (qsat-RS) obtained from RS pressure and temperature profiles, and the saturation specific humidity profile (int-q-RS) obtained by integrating RS vertical gradient of refractivity (M_RS) are shown. Panels (a), (b) and (c) are for the high mode (vertical resolution 375 m, interpolated every 150 m), while panels (d), (e) and (f) are for the low mode (vertical resolution 150 m).

Another improvement consisted in testing different values of the Hlim level in case of the presence of relative maxima in radar reflectivity profiles. An example is provided in Fig. 2. The $C_n^{\mathrm{2}}$ profile in Fig. 2a shows three peaks. The dominant peak at 902 m corresponds to the lower level of a cloud, since q observed with the RS is saturated (the red solid and dashed lines are superimposed in Fig. 2b). The double integration of the radar data upward and downward to this level leads to the radar profile of q illustrated in Fig. 2b (black line). Radar and RS profiles agree up to 1500 m, but above this level the radar and RS profiles deviate. In this specific case the downward integration fails (below 3100 m). We recall that in the upper part of the profile, the integration is performed from 4800 m down to Hlim = 902 m. The combination of $C_n^{\mathrm{2}}$, ϵ and shear measured by the radar fails to reproduce the discontinuity present in the RS profile between 3100 and 2800 m. This may be due to the spatial heterogeneity of the air mass, but there is no specific ancillary information to prove it. Different values of Hlim (1802 and 3452 m) were also tested, based on the consideration of different peaks of the $C_n^{\mathrm{2}}$ profile. The best result is obtained when Hlim = 1802 m, which is the height of the second peak from ground. The result is illustrated in Fig. 2c. It shows a slight improvement of the agreement in comparison to results shown in Fig. 2b, although the divergence of the radar profile below 3100 m down to 2400 m is still obvious. The transition at 1802 m does not correspond to any marked transition in the RS humidity profile. This choice probably improved the result due to its central location in the profile.

Figure 2d to f illustrate the results for the same profile, when the low mode is used instead of the high mode. In this case, the vertical resolution of the radar measurements is 150 m instead of 375 m. This can be checked on the $C_n^{\mathrm{2}}$ profile in Fig. 2d, which is not as smooth as in Fig. 2a. The agreement between the radar and RS profiles of q is slightly improved, but again not striking (Fig. 2f compared with Fig. 2c).

Figure 3 Specific humidity profiles obtained for two different profiles from the BLLAST data set, using an automatic detection of the initial specific humidity value (q_o) at the level z_o, where the integration is initialized (a, c) and through an adjustment of q_o (b, d). The vertical resolution is 75 m. Same details as in Fig. 1c.

Finally, another issue was raised by the choice of the initial conditions. We constrained the value of q_o to be equal to the RS observation at the same level. This level sometimes coincided with a level of sharp humidity decrease, as illustrated in Fig. 3a at 3000 m or in Fig. 3c at 4400 m. This was frequently the case at the upper boundary during BLLAST since the range gate where the radar echoes vanished also coincided with a sharp decrease of moisture. So we rectified the initial value by hand, to avoid an accumulation of the error along the whole profile, as shown in Fig. 3, where the humidity profiles prior (left column) and posterior (right column) to the correction are compared. We are well aware that this manual correction could not be applied for an operational purpose, but attempts to introduce this correction procedure in a dedicated algorithm are currently underway.

Figure 4 (a, c, e, g) Scatterplots of radar versus RS specific humidity during June and July 2011 (BLLAST), September 2012 (HyMeX SOP1), October 2012 (HyMeX SOP1) and February 2013 (HyMeX SOP2), respectively, with the linear regression line (in red) and the 1:1 slope line (in black). The R² correlation coefficient of the regression is also specified. (b, d, f, h) Deviation profiles between the RS specific humidity measurements and radar-based estimate for the above specified data sets ± the standard deviation. The mean bias and standard deviation values for the whole dataset are reported, along with the maximum standard deviation value (g kg⁻¹).

4.2 Comparison between the RS and radar retrieved humidity profiles

We gathered in Fig. 4 the results of the comparison between the radar estimation of the humidity profile and the observed RS humidity profiles for the three data sets in low mode. Panels (a), (c), (e) and (g) show the scatter plots of radar vs. radiosonde data, and panels (b), (d), (f) and (h) show the vertical profile of the deviations between radiosonde and radar data, all of them obtained from low mode data. The red line in the scatter plots represents the linear regression line drawn from the data, which can be compared to the slope 1:1, i.e., the black line. HyMeX SOP1 data, the most numerous in terms of number of profiles, were split into two graphs. The vertical resolution of the BLLAST WPR measurements is 3 times higher than the one characterizing HyMeX measurements, which explains the larger sample size in Fig. 4a. Results are also summarized in Table 1, where the high mode has been added.

The best results are clearly those obtained during HyMeX SOP2, since the correlation coefficient is exceeding 0.93, the mean bias is 0.04 g kg⁻¹ or smaller and the mean standard deviation is 0.18 g kg⁻¹ or smaller. These successful retrievals are linked to the peculiar conditions of the HyMeX SOP2 experiment (limited variability of the humidity content). So we cannot rely upon this example to assess the effectiveness of our method. In most other cases, with the exception of the BLLAST data in high mode, the correlation coefficient exceeds 0.8, the slope of the regression lines is close to unity, the mean bias is equal or smaller than 0.24 g kg⁻¹ and the standard deviation is not exceeding 1 g kg⁻¹ (Table 1).

Table 1 Correlation coefficient R² of the linear regression between radar and RS humidity values for the three data sets. Mean bias (RS minus radar) and standard deviation for the whole dataset are also specified, together with the largest standard deviation (SD) value in the dataset. “lm” and “hm” stand for low mode and high mode, respectively.

We identified the reason why the BLLAST dataset gave poorer results in high mode than in low mode. In high mode these are 0.73, 0.36 g kg⁻¹ and 1.19 g kg⁻¹ for the correlation coefficient, mean bias and standard deviation, respectively (Table 1). These relatively lower-quality results are due to the presence of thin layers of dry air (−300 m vertical depth) that the radar fails to reproduce in high mode since its vertical resolution (375 m) is not high enough.

The other results in Table 1 do not exhibit an outperformance at any of the two modes. Concerning the variability of the error and its standard deviation along the vertical, we cannot draw any general conclusion even if the HyMeX results seem to be better (smaller bias, smaller standard deviation) above 2500 m (Fig. 4d and f). The bias remains small and similar for BLLAST within the whole profile.

To conclude, we consider that the method we propose yields good radar profiles to start the processing at a finer time resolution. Before getting the intermediate profiles, let us examine the variety of calibration coefficients α² that were obtained for the three data sets.

4.3 Variability of the calibration coefficients

With the aim of comparing our estimates of the calibration coefficient α² in Eq. (10) with values found in literature and studying its variability relative to the stability conditions, we first considered the coefficients within each dataset, which enabled us to avoid possible issues about the radar calibration. According to Eqs. (8) and (10) and considering the same WPR, the variability in α² should reflect the variability of the ratio of the two outer scales $\frac{L_w}{L_{\mathit{\varphi }}}$. Gossard et al. (1998) expected this quantity to be dependent on the stability conditions, at least inside the ABL, where the turbulence is homogeneous. Small ratios should occur under stable conditions, which would correspond to large α² values, and the opposite behavior under unstable conditions. To check this hypothesis, we used the radiosonde data to get an estimation of the stability conditions based on the gradient Richardson number:

$$\begin{array}{}\text{(19)}& R_i={\displaystyle \frac{g}{{\mathit{\theta }}_v}}{\displaystyle \frac{\frac{\mathrm{d}{\mathit{\theta }}_{\mathrm{v}}}{\mathrm{d}z}}{{\left(\frac{\mathrm{d}V}{\mathrm{d}z}\right)}^{\mathrm{2}}}},\end{array}$$

where θ_v is the virtual potential temperature. We chose the BLLAST dataset, since summer conditions provide the largest range of stability conditions. Figure 5a shows the calibration coefficients obtained from the BLLAST dataset in low mode as a function of the stratification with a differentiation of the lower part and upper part of the profiles, located below and above Hlim, respectively. The color scale indicates the hour of the day. During the (dry) convective period of the day, the lowest portion of the profiles is expected to correspond to a mixed ABL, with Richardson numbers potentially negative in the surface layer and close to zero above. The corresponding calibration coefficients can be clearly identified in the bottom left corner of Fig. 5a. Smaller values (below 0.01) are observed under unstable conditions, in accordance with Gossard et al. (1998). However, a significant number of the coefficients corresponding to unstable conditions are also located in the range of 0.03–0.5, where most of the coefficients of the lower layer can be found (the logarithmic mean of α² in the lower layer is 0.11). In contrast, nighttime coefficients measured within the lower layer (dark blue and orange circles) are preferentially large (> 0.5), while the Richardson number can vary but is usually close to or larger than 0.25, i.e., the critical value of the Richardson number, indicating stable conditions. As expected, the upper layers are slightly or clearly stable (squares). In this case the logarithmic mean is 0.16, not far from the value observed in the lower layer. According to VanZandt et al. (1978), the variability of the coefficients in the upper part of the profiles is linked to the variation of the filling factor and depends on the lapse rate of the free troposphere. Anyhow, this variability (the squares span roughly 2 orders of magnitude) is less marked than the variability due to the stability conditions of the low layers (the circles span 4 orders of magnitude). The average value of the calibration coefficients for BLLAST is 0.13, which is close to the lower boundary of coefficients [0.26–1.11] proposed by Gossard et al. (1998).

Figure 5 (a) α² calibration coefficients obtained with the BLLAST dataset in low mode as a function of the stratification estimated with the gradient Richardson number. The results are divided into the lower (circles) and upper part (squares) of the ABL and shown as a function of the time of the day (colored scale). (b) Scatterplot of the calibration coefficients α² (high level versus low level) obtained with the BLLAST and HyMeX datasets in low mode. The crosses correspond to the profiles measured during the potentially dry convective period, i.e., 09:30–17:00 UTC.

The BLLAST coefficients obtained with the high mode are not reported since we estimated they are not representative, especially in the lower part of the profiles. First, the average of the transition levels (Hlim) is 1193 m a.g.l. in high mode versus 928 m a.g.l. in low mode, and the distribution is clearly shifted towards higher Hlim values in high mode, which means that the transition level is certainly higher than the boundary layer top in most cases. Second, we obtained rather large values of the coefficients (6 values larger than 10) for the lower part of the profiles, which we attribute to the poorer vertical resolution in high mode in an area of sharp variations of the humidity.

The same analysis was applied to the HyMeX datasets from which no clear result arose, during either SOP1 or SOP2 (not shown). During HyMeX, the development of the boundary layer was most of the time generated by mechanical turbulence (due to the wind intensity or to the roughness change at the sea–land transition). This can explain why the gradient Richardson number is not a good indicator of the variability in the calibration coefficients under the HyMeX conditions. There is also no clear difference in the HyMeX calibration coefficient values between the lower and upper parts of the profiles, probably because moist convection equally affected all levels. Although the HyMeX radar was calibrated (in low mode), the coefficient values determined during SOP1 were 3 times smaller than those found during BLLAST. We artificially shifted the HyMeX calibration coefficients (by multiplying them by 3) to make the logarithmic average match in order to be able to compare the variability. Results are illustrated in Fig. 5b in terms of upper layer versus lower layer α² coefficients. As seen before, due to the variability in stability conditions, BLLAST is characterized by a larger span of lower layer coefficient values. The calibration coefficient values determined during HyMeX SOP1 and SOP2 span roughly 2 orders of magnitude in the lower layers and 2 orders of magnitude in the upper layers on each side of the 1:1 line.

The main conclusion we can draw from these results, coming primarily from the BLLAST data, is the necessity of distinguishing between the mixed layer and the free troposphere in case of unstable conditions in the low troposphere. Figure 5b also suggests that the filling factor proposed by VanZandt et al. (1978) for the upper layers does not vary that much. This explains why some authors obtained satisfying results by using a constant calibration coefficient value with VHF-band WPRs when sensing this portion of the atmosphere.

We propose an explanation for the different calibration coefficient values obtained during HyMeX SOP1 (0.04 on average) and BLLAST (0.13 on average). We recall that both radars were calibrated. We assume that we can rely on this calibration since (i) values of $C_n^{\mathrm{2}}$ determined for BLLAST in high mode are similar to those obtained in low mode and (ii) the two calibrations, based on the rain gauge measurements, are independent. During HyMeX SOP1, atmospheric thermodynamic conditions were close to saturation most of the time (as shown for example in Fig. 2c). Although we discarded the rainy radar profiles when the rainfall affected the whole air column, it is likely that isolated pockets of rainfall may have locally increased the radar $C_n^{\mathrm{2}}$, which implies a decrease in α². We checked the distributions of $C_n^{\mathrm{2}}$ for the three datasets and found that the logarithmic averages of $C_n^{\mathrm{2}}$ (close to the median values) are 1.4, 31 and 1.0 × 10⁻¹⁴ m${}^{-\mathrm{2}/\mathrm{3}}$ for BLLAST, HyMeX SOP1 and HyMeX SOP2, respectively, which confirms the former hypothesis of higher $C_n^{\mathrm{2}}$ values (and consequently lower calibration coefficient values) for the moist conditions during HyMeX SOP1. This disturbance can in principle be considered as a limit in the application of the method we propose, since the clear-air turbulence conditions required for the application of Eq. (7) is not totally fulfilled. However, we consider that, based on the good results we obtained (see Fig. 4c and e), the method can be successfully applied to HyMeX SOP1.

Table 2 Parameters used to retrieve 15 min spaced vertical profiles of q in the time interval between two radiosoundings RS1 and RS2. Hlim serves to apply the convenient calibration coefficient and to merge the two integrals.

5 Continuous humidity monitoring between radiosonde observations

5.1 Method

Our objective in this section of the paper is to check whether radar data can be successfully used to describe the detailed structure of the low troposphere between two consecutive radiosoundings, even if the latter are 6 or 12 hourly. The two bordering soundings, called RS1 and RS2, are used as initial and final conditions. We first retrieved the radar humidity profiles closest to RS1 and RS2 on the basis of the application of the algorithm described in the previous sections. The algorithm provides the values of Hlim and the calibration coefficients at the launching times of RS1 and RS2. RS1 and RS2 also provided the initial and final bottom and top border conditions for q, which we interpolated at the radar times (every 15 min). Pressure and temperature values were also interpolated at the times between RS1 and RS2 (at all levels) to provide θ and P present in Eq. (6) and also to constrain q to the humidity saturated value based on the interpolation of the RS data at the same time and level. In addition, the sign of M at each level was taken as the sign of the RS1 humidity gradient at the same level for the first half of the period separating RS1 and RS2 and as that of RS2 for the second part of the period. Finally, the radar calibration coefficients α² were time-interpolated. In contrast, values of Hlim obtained from the radar data at the launching times of RS1 and RS2 were not interpolated. The radar provided an updated value of Hlim and an updated vertical profile of M² every 15 min (through new $C_n^{\mathrm{2}}$, shear stress and ϵ profiles). Ultimately, using Eq. (6), the vertical profile of q may be retrieved at a fine time resolution of 15 min. The choices for the parameters just described are summarized in Table 2.

In case of several values of Hlim due to the presence of several peaks in the $C_n^{\mathrm{2}}$ profile obtained at time t, a continuity criterion was used to select the appropriate value. This criterion was applied to a variety of HyMeX case studies, which frequently showed multiple layers with sharp humidity gradients. In contrast, the $C_n^{\mathrm{2}}$ profiles observed during BLLAST exhibited a marked isolated peak that could be directly taken as Hlim.

Figure 6 Vertical profiles of humidity measured during HyMeX on 24 September 2012: the radar and RS profiles used as time border conditions are shown in panels (a) and (b) with the same details as in Fig. 1. These border conditions are reproduced in panels (c) and (d) with the thick blue solid lines for RS1 (03:14 UTC) and red one for RS2 (09:00 UTC). Intermediate radar profiles (thin solid lines) and corresponding saturated water vapor profiles (dashed lines) are presented in panel (c). Panel (d) shows the lidar humidity profiles (thin lines with crosses), the radar ones (thin lines) and the humidity profiles resulting from a linear interpolation between RS1 and RS2 (dashed lines).

5.2 Some results

A first example of the results obtained during HyMeX is shown in Fig. 6. RS1 and RS2 are 6 hourly and were chosen to illustrate a case study, on 24 September 2012, when the moist lower troposphere dries and a mixed boundary layer develops. The RS1 humidity profile at 03:14 UTC is close to saturation from 500 up to 2000 m (Fig. 6a), whereas the RS2 profile at 09:00 UTC is drier and shows a mixed layer with a depth of 1350 m and a constant q value of 8.4 g kg⁻¹ (Fig. 6b). The humidity profiles retrieved by the radar between 03:14 and 09:00 UTC are represented with thin solid lines in Fig. 6c, with the time being color coded. The dashed lines are the corresponding saturated profiles calculated from the profiles of P and T obtained by interpolating the data from the two RS. The radar profiles gradually dry in the lower levels and a mixed boundary layer develops from 07:14 to 09:00 UTC, accompanied by a decrease in air temperature (with a consequent decrease of q saturated between RS1 and RS2 in Fig. 6c). In fact, the weak low-level wind, which had blown from the south in the early morning, turned gradually to northwesterly wind between the sunrise (around 05:00 UTC) and 06:14 UTC (not shown). The mist marine layer was then replaced by a continental, cloud-free boundary layer that could easily develop due to an increase in the wind strength (15–20 m s⁻¹) and to a larger surface–air temperature contrast. The radar was particularly helpful to detect the top of the mixed layer.

In Fig. 6d, we compare the radar humidity profiles (solid lines) to the 15 min averaged profiles obtained with the lidar at the same times (crosses) and to the humidity profiles calculated from a linear interpolation of the RS1 and RS2 profiles (dashed lines). The first four lidar observations (between 03:14 and 04:59 UTC) are attenuated above 750 m because of the presence of a cloud layer. In fact, atmospheric particles can lead to antagonistic effects on the lidar beam: few and scattered particles may lead to an increase of the backscattered radiation, whereas dense particle ensembles, as those found in a thick cloud, are usually characterized by large optical thicknesses (> 1–2), which translates into laser beam attenuation overwhelming particle backscattering. In the early morning, the sharp decrease of the lidar signal above 750 m in Fig. 6d clearly reveals the base of a thick cloud at 750 m. As the air dries up with the time the lidar recovers its capacity to cover the lower troposphere and the final lidar humidity profiles are very close to the RS2 and radar profiles.

Figure 7 Time–height cross sections of q observed on 24 September 2012 by the lidar (a), calculated from the radar data (b) and interpolated between RS1 and RS2 (c). The dots are for the radar saturated values in panel (b) and RS saturated values in panels (a) and (c).

The representation in Fig. 6d hardly enables us to distinguish between the radar humidity and the humidity calculated from the RS interpolation. That is why we preferred to show the same results with height–time cross sections of humidity (Fig. 7). In this figure, we highlight with black dots the levels characterized by saturation conditions. Due to the difficulty to rely on the lidar humidity profiles measured under saturated conditions, we superimposed on the lidar map in Fig. 7a, the dots obtained with the RS data (the same as those in Fig. 7c). The thick cloud remains over the lidar until 05:00 UTC. Simultaneously, the radar detects saturated values of q (Fig. 7b). After a short period of clear air, another cloud is advected over the measurement site between 05:30 and 06:30 UTC, with a base at 1100 m, as indicated by the lidar. This cloud is also captured by the radar, with saturated values of q from 1500 to 2200 m, at 05:30 and 06:30 UTC. After 05:00 UTC the lidar shows humidity values similar to those measured by the radar and the RS in the low layers, although the mixing of the boundary layer is best represented by the radar and probably the most likely since it results from an accurate measurement of the inversion height, which is well marked by the Hlim transition. The moisture retrieved by the radar exhibits a mixed boundary layer, whose depth increases from 500 (at 07:00 UTC) to 1300 m (at 09:00 UTC) and dries with the time.

Figure 8 Radar-based humidity profiles and humidity height–time cross sections for the 23 September HyMeX case study: from RS1 (09:02 UTC, blue solid line in a) to RS2 (14:59 UTC, red solid line in a), from RS2 (blue, b) to RS3 (20:58 UTC, red solid line in b) and from RS1 (blue, c) to RS3 (red, c). First row is for the radar profiles (same details as in Fig. 6), second for the height–time humidity of the lidar, third for the radar and fourth for the interpolated RS. The dots in the last three rows demarcate the saturated values as in Fig. 7.

With the following example (case study on 23 September 2012 from HyMeX SOP1), we extended the time separating the two border radiosounding by choosing RS1 at 09:01 UTC and RS3 at 20:31 UTC. This enabled us to use an intermediate RS at 14:58 UTC, called RS2, and to compute the radar profiles during three distinct periods, namely 09:01–14:58 UTC (Fig. 8a, d, g and j), 14:58–20:31 UTC (Fig. 8b, e, h and k) and 09:01–20:31 UTC (Fig. 8c, f, i and l). The juxtaposition of the first two columns should give the third. As expected, this is clearly the case for the lidar (Fig. 8d, e and f) but is not so obvious for the interpolated RS (Fig. 8j, k and l).

The RS bottom border conditions in Fig. 8l do not recreate the drying that occurs at 500 m between 09:00 and 15:00 UTC (Fig. 8a and j; 12 down to 7 g kg⁻¹) and the following moistening from 15:00 to 20:30 UTC. The saturated areas (black dots) in panels (l) or (j) and (k) are not consistent either. Remember that dots in this row indicate the interpolated RS. In fact, the decrease in q at 500 m in panel (j) widens the gap to the saturation, especially as the saturation values climb due to the increase of temperature with time. In contrast, the low-layer moistening from 15:00 to 20:30 UTC rapidly leads to saturation conditions, not at the surface where the air is warm, but at a higher altitude of 1300 m (Fig. 8k around 16:00 UTC). A deeper cloud layer appears around 18:00 UTC as testified by the RS interpolation in panel (k). The thickening of the cloud layer occurs later, around 19:30 UTC in panel (l). Consequently, a simple interpolation within the 12 h interval between RS1 and RS3 seems unrealistic since it does not reflect the proper daytime evolution of the low-level boundary conditions. This also influences the saturation conditions at the mid-level.

Figure 9 Same as in Fig. 7 for 26 September 2012.

The comparison with the lidar enhances the difficulty for the radiosoundings to describe the humidity fields under non-steady atmospheric conditions. On 23 September 2012, a strong convective line was active, far west of the site, extending from Iceland to the west of the Iberian Peninsula (this convective line finally crossed the measurement site during the following night). This frontal activity generated a pocket of moist air in the lower layers that was advected by southeasterly winds (150^∘ at 500 m) over the western Mediterranean sea, while it circumvented the Pyrenean mountains. Another pocket of moist air was situated to the east of the site due to some convective activity over Italy. Between these two moist areas, there was a pocket of drier air that moved according to the relative influence of the two convective areas on each side. Consequently, the measurement site, which was located in this area, encountered varying moisture conditions.

The lidar captures well the large humidity variability characterizing the lower levels with large amounts of q (13–14 g kg⁻¹) in the time interval 09:30–13:30 UTC (Fig. 8d), as a result of the influence of the western convective line, and the drying of the lower layers in the time interval 13:30–20:30 UTC (Fig. 8e), when the moisture pocket moved further to the west. The 09:00 UTC RS that was launched during a period of increasing moisture (12 g kg⁻¹ between 500 and 800 m) is not able to capture the whole increase (unfortunately the lidar data were missing at 09:00 UTC as indicated by the blue vertical stripe from 09:00 to 09:45 UTC in Fig. 8d). Consequently the radar, whose initial conditions are based on the 09:00 UTC RS, fails to capture the large humidity amounts characterizing the lower layers between 09:30 and 13:30 UTC (Fig. 8g).

Between 13:30 and 17:00 UTC, the radar estimations below 1500 m are closer to the lidar measurements. Both instruments indicate a top level of the moist layer varying between 1200 and 1500 m (Fig. 8d–e and g–h). Above this level, the lidar beam is extinguished by the thick clouds. The detection of this transition level is facilitated, for the radar, by its being marked by a large wind shear which favors the $C_n^{\mathrm{2}}$ increase.

Between 17:00 and 18:30 UTC, the radar and the lidar both detect higher humidity amounts in an intermediate layer between 800 and 1500 m or 500 and 1100 m, respectively (Fig. 8e and h). Within this layer, the radar shows saturated values close to the levels where the interpolated RS values are saturated, but limited between 17:00 and 18:30 UTC. The radar humidity values are also consistent with the lidar values. Between 16:00 and 18:30 UTC, the radar detects a cloud in the vertical region 2000–2500 m (Fig. 8b and h). A cloud is also well visible in the particle backscatter field (not shown) obtained by the lidar. The particle backscatter data at 1064 nm are able to properly reveal both aerosol layers and cloud/precipitation particles. Specifically, hydrometeors evaporating/sublimating before reaching the ground are observed (as vertical thin stripes) between 17:30 and 19:00 UTC in the vertical region 0.5–1.2 km.

Finally, if we now consider the capacity of the radar to retrieve the humidity variability during the whole period from 09:00 to 20:30 UTC (Fig. 8i), we must recognize that the radar fails to reproduce the lidar data (Fig. 8f), with the exception of the transition level around 1400 m, between 13:00 and 20:30 UTC. The cloud mentioned previously is unfortunately not seen, probably due to the erroneous calibration coefficients imposed by the too large interval separating the two bordering RS.

The 23 September 2012 case study showed that, for changing conditions near the surface (that are not typical of a classical diurnal cycle evolution) and 12-hourly radiosoundings, both the combined radar–RS algorithm and the simple time interpolation between two radiosoundings fail to reproduce the detailed structure of the humidity in the low troposphere, essentially because of a lack of documentation of the evolution of the humidity at the first radar gate. As stated before, the data from a ground station would not have been helpful due to the variability between the surface and the first radar gate (see for instance Fig. 8e). In contrast, when the time interval between the two RS was reduced to 6 h, the radar proved to be able to retrieve the humidity amounts measured by the lidar and the presence of a cloud in the intermediate layers.

Another example to illustrate the performance of the radar-based approach with respect to a RS interpolation and the lidar measurements is shown in Fig. 9. RS1 and RS2 on 26 September 2012 are 6 hourly. Both show clear air and a regular reduction of the humidity content with the height (Fig. 9c). In fact, a pocket of moister air occurred between the two RS at 03:00 and at 09:00 UTC that is well observed with the lidar (Fig. 9a) and partly retrieved by the radar (Fig. 9b), with local saturated values. Even if the radar retrieval shows a few erroneous profiles (04:30, 05:30, perhaps 07:00 and 07:30 UTC, etc.), the radar documentation is necessary to detect the morning moistening of the air.

6 Conclusions

We demonstrated in the first four sections of this paper that although WPRs, with their first three moments, measure essential parameters for the determination of the vertical humidity gradient, the radar data cannot be used to retrieve the vertical profiles of humidity independently from other sensors' data. To obtain the profiles, we applied a method already proposed by Tsuda et al. (2001), which consists in using a combined retrieval algorithm exploiting WPR measurements supported by RS observations at a coarser frequency. This algorithm is based on several approximations and assumptions that proved to be appropriate since the accuracy of the results we obtained did not exceed 0.25 g kg⁻¹ (mean bias between q radar and q RS). To obtain these results, we improved the algorithm proposed by Tsuda et al. (2001) by using a key parameter from the radar, which is the level of the reflectivity peak value, Hlim, to split the calculations in two parts, with two different calibration coefficients accounting for two distinct vertical regions of different turbulence characteristics. The introduction of this level also mitigated the errors by replacing a long integration by two shorter ones.

After assessing the algorithm at the time of the RS observations, we applied it between two RS profiles to obtain humidity profiles at a finer time resolution and to check the performance of the combined algorithm with respect to a simple RS time interpolation. We used, when available, simultaneous lidar data to assess the results. The set of data that enabled this comparison was collected during a period seldom characterized by the presence of clear-sky conditions, while cloudy conditions prevailed (HyMeX SOP1). In the presence of clouds, the lidar beam is rapidly attenuated above cloud base, so that the assessment can only be made in the lower portion of the profiles.

We obtained some satisfactory results, provided that the time separating the two boundary RS did not exceed 12 h. However, we also encountered some hindrances that make the method hard to apply in an entirely automatic way due to the assumptions we made. These difficulties are summarized below:

• The most restrictive issue is the one associated with the border conditions (bottom and top). The method assumes that they vary linearly between the two RSs, which is not always true. If the border conditions are not well defined (for instance at a level of strong moisture gradient), the error may propagate and become large at the Hlim level. Additionally, the resulting profiles can easily move apart, towards the two constraining borders: either towards 0 g kg⁻¹ as the minimum value or towards the saturated moisture content as the maximum value.

• Although pressure and temperature are secondary parameters in the algorithm, so that their estimate does not need to be as accurate as the border conditions, the profiles for these parameters have to be provided. These two parameters are also used to constrain the computed humidity values to the saturation q. We used a linear interpolation of the two border RS to get the intermediate P and T profiles. Alternatively, these profiles could be provided by models, which are usually more reliable for pressure and temperature than they are for humidity.

• The constraint on the sign of the humidity gradient is also an issue that can hardly be solved by a simple interpolation or a continuity constraint in time. Some authors constrain their results with GPS measurements of the integrated water content. This approach failed with our data set.

• We highlighted the necessity of calibrating the vertical gradient of radar refractivity, with calibration coefficients likely to vary in time and space. This revealed to be helpful, but in this case a simple interpolation between the initial and final coefficients could be too large an approximation. However, the detection of the transition level between the boundary layer and the free troposphere was definitely helpful.

Finally, we demonstrated that the combined RS–radar algorithm used to retrieve the humidity profiles outperforms a simple interpolation of the RS observations. The radar is especially skilled at determining the evolution of the transition layers, which is usually an issue when using other remote-sensing measurements such as radiometer measurements. However, the present method should be used with caution and is probably more adequate for post-processing of a dataset for scientific purpose than for a blind use in an automatic platform.