4-D-VAR assimilation of disdrometer data and radar spectral reflectivities for raindrop size distribution and vertical wind retrievals

Abstract. This paper presents a novel framework for retrieving the vertical raindrop size distribution (DSD) and vertical wind profiles during light rain events. This is also intended as a tool to better characterize rainfall microphysical processes. It consists in coupling K band Doppler spectra and ground disdrometer measurements (raindrop fluxes) in a 2-D numerical model propagating the DSD from the clouds to the ground level. The coupling is done via a 4-D-VAR data assimilation algorithm. As a first step, in this paper, the dynamical model and the geometry of the problem are quite simple. They do not allow the complexity implied by all rain microphysical processes to be encompassed (evaporation, coalescence breakup and horizontal air motion are not taken into account). In the end, the model is limited to the fall of droplets under gravity, modulated by the effects of vertical winds. The framework is thus illustrated with light, stratiform rain events. We firstly use simulated data sets (data assimilation twin experiment) to show that the algorithm is able to retrieve the DSD profiles and vertical winds. It also demonstrates the ability of the algorithm to deal with the atmospheric turbulence (broadening of the Doppler spectra) and the instrumental noise. The method is then applied to a real case study which was conducted in the southwest of France during the autumn 2013. The data set collected during a long, quiet event (6 h duration, rain rate between 2 and 7 mm h−1) comes from an optical disdrometer and a 24 GHz vertically pointing Doppler radar. We show that the algorithm is able to reproduce the observations and retrieve realistic DSD and vertical wind profiles, when compared to what could be expected for such a rain event. A goal for this study is to apply it to extended data sets for a validation with independent data, which could not be done with our limited 2013 data. Other data sets would also help to parameterize more processes needed in the model (evaporation, coalescence/breakup) to apply the algorithm to convective rain and to evaluate the adequacy of the model's parameterization.

Abstract. This paper presents a novel framework for retrieving the vertical raindrop size distribution (DSD) and vertical wind profiles during light rain events. This is also intended as a tool to better characterize rainfall microphysical processes. It consists in coupling K band Doppler spectra and ground disdrometer measurements (raindrop fluxes) in a 2-D numerical model propagating the DSD from the clouds to the ground level. The coupling is done via a 4-D-VAR data assimilation algorithm. As a first step, in this paper, the dynamical model and the geometry of the problem are quite simple. They do not allow the complexity implied by all rain microphysical processes to be encompassed (evaporation, coalescence breakup and horizontal air motion are not taken into account). In the end, the model is limited to the fall of droplets under gravity, modulated by the effects of vertical winds. The framework is thus illustrated with light, stratiform rain events.
We firstly use simulated data sets (data assimilation twin experiment) to show that the algorithm is able to retrieve the DSD profiles and vertical winds. It also demonstrates the ability of the algorithm to deal with the atmospheric turbulence (broadening of the Doppler spectra) and the instrumental noise. The method is then applied to a real case study which was conducted in the southwest of France during the autumn 2013. The data set collected during a long, quiet event (6 h duration, rain rate between 2 and 7 mm h −1 ) comes from an optical disdrometer and a 24 GHz vertically pointing Doppler radar. We show that the algorithm is able to reproduce the observations and retrieve realistic DSD and vertical wind profiles, when compared to what could be expected for such a rain event.
A goal for this study is to apply it to extended data sets for a validation with independent data, which could not be done with our limited 2013 data. Other data sets would also help to parameterize more processes needed in the model (evaporation, coalescence/breakup) to apply the algorithm to convective rain and to evaluate the adequacy of the model's parameterization.

Introduction
Knowledge of the raindrop size distribution (DSD), both at the ground and throughout the vertical atmospheric profile, is an important topic. In remote sensing, parameters of the Z-R relationship, allowing a radar reflectivity to be converted into a rain rate, depend on rain microphysics. These parameters are generally based on a particular DSD distribution and are supposed to be constant through the entire atmospheric column. In practice, rain rates estimated from reflectivity through the Z-R relationship can vary at least by a factor of 2 depending on the DSD (List, 1988). In telecommunications, and especially for ground to space telecommunication, rain-induced attenuation also depends on rain microphysics, especially in the Ku band and beyond. More generally the observation of vertical DSD profiles is fundamental to address many questions regarding the numerous processes involved during a rain event. A large number of investigations have been conducted in order to model the vertical Published by Copernicus Publications on behalf of the European Geosciences Union. evolution of the DSD during rain for different meteorological situations (List and McFarquhar, 1990;Brown Jr, 1993;Hu, 1995;Prat and Barros, 2007a;Prat and Barros, 2007b;Barros et al., 2008;Mcfarquhar, 2010). The parameterization of the different phenomena occurring during droplet fall is a complex task due to the great number of processes involved. Among them, some are considered as predominant, notably a raindrop can break up, coalesce or evaporate. It is also subject to sorting by gravity along with updrafts, downdrafts and horizontal wind. All these processes, playing an important role in the nature of the DSD, are the subject of numerous studies. More specifically, raindrop sorting induced by vertical drafts can play an important role in the vertical DSD variability. For instance, Kollias et al. (2001) show that updraft structures can cause horizontal and vertical sorting of raindrops with a lack of large raindrops (> 3 mm) in the updraft core and an increase at the periphery. Moreover, concerning evaporation, drizzle can never reach the ground in specific meteorological situations (VanZanten et al., 2005). On the other hand, Mcfarquhar (2010) summarizes the main results on models concerning the coalescence and breakup processes. However, he explains that these models, built on laboratory or simulation experiments, generally suffer from a lack of validation under real conditions. Nevertheless, a very large number of direct or indirect observations of DSD are available either on the ground, from disdrometers (among others, Joss and Waldvogel, 1967;Kruger and Krajewski, 2002;Delahaye et al., 2006), or at various heights, from the use of Doppler radars (Kollias et al., 2002;Peters et al., 2005). Several studies have been carried out to compare these different observations (Bringi et al., 2003;Tokay et al., 2009;Schleiss and Smith, 2015). However, the combination of these sources of information about the DSD is a difficult task given their very different natures. It is even more demanding when you add the different spatiotemporal scales and the various locations of the measurements.
In order to improve vertical DSD profile retrievals, this study establishes a framework allowing the assimilation at a fine scale (< 5 min, 100 m) of heterogeneous observations in a dynamic model. The developed method allows the use of observations from different instruments, each with its own spatiotemporal resolution. Compared to direct observation retrievals, merging data in a dynamic model includes spatiotemporal consistency to the DSD retrievals. Moreover, it could help to assess and improve the model parameterization.
In this paper, we will focus on disdrometer and vertical radar reflectivity collocated observations. The data are merged through the use of a 4-D-VAR assimilation algorithm. The evolution in time and space of the DSD is based on a very simplified version of the Hu and Srivastava (1995) propagation model (see Eq. 3). It is actually limited to the fall of droplets under gravity, solely modulated by the effect of vertical wind. The goal of this study is to develop a new mathematical framework in a simplified context. Once introduced, it will be tested by evaluating its performance, sensitivity and limits on both simulated data and real data (a calm stratiform real event of light rain). Note that, even if it is simple, this model ensures a DSD which matches the measurements of all instruments (and their different scales). This aspect is not negligible when, as mentioned previously, you consider all the work done to establish consistent Z-R relationships.
Improvements of this framework (encompassing evaporation, horizontal air motion, coalescence/breakup) as well as validation, both requiring more measurements, are deferred to future work allowed by extended data sets.
The paper is organized as follows. In Sect. 2, we present the simple model used and discuss the terms of the complete Hu and Srivastava (1995) propagation model, namely wind, collision and evaporation terms. Then, model simplifications are retained considering different hypotheses. General principles behind 4-D-VAR data assimilation are presented and a focus is made on the cost function and regularization terms. Disdrometer and radar observation operators are then given. Finally we describe the algorithm used to retrieve the DSD and vertical wind fields from radar and disdrometer data. Then, in Sect. 3, the instruments available and the case study are described. In Sect. 4, through a twin experiment (i.e., on simulated data), we explore the impact of observation errors in addition to non-modeled phenomena on the algorithm performances as well as the quality of the retrievals. Section 5 gives some results obtained in a stratiform real case study. Finally, Sect. 6 concludes the study. Perspectives are also drawn concerning the improvements needed to extend the application of the method to a more general context.

Assimilation method
The aim of this work consists in retrieving information about the drop size distribution (DSD) vertical profiles by linking the measurements made by different instruments at different scales and heights. To carry out this relationship, we will use a dynamic model which propagates the DSD through space and time (denoted below as "propagation model"). In this section, firstly, we describe the corresponding partial differential equation (PDE) model used. We also introduce the associated discretization scheme, and the unknowns of our model. Then, we explain how the 4-D-VAR data assimilation algorithm combines the data available through the model to retrieve its unknowns. The second part of this section details this algorithm.

Propagation model
As explained before, this study develops a simple framework and so uses a very simple model presented in Sect. 2.1.1. This simple model has the merit to add a spatiotemporal coherence to the DSD field. Since it only incorporates the effects of gravity and vertical wind, it lacks the complexity required to fully model rain. In Sect. 2.1.2, we introduce a complete PDE governing the evolution of the DSD (Hu and Srivastava, 1995). This propagation model is not within the purview of this article. It is solely used to discuss the different terms of the complete PDE (Eq. 3) and see how their inclusion or exclusion could affect the results. Finally, in Sect. 2.1.3, we present the effective "simplistic" numerical model and the discretization underlying hypotheses.

Simple propagation model
The DSD, noted N, is the number of raindrops by unit of volume and diameter (unit: m −4 ). It is a function of time, position and diameter. In the framework presented here, we will work in an air column, and so limit the study to a timeheight-diameter space, so that N = N (t, z, D). Then, the PDE modeling the vertical evolution of DSD has the form: The first term of Eq. (1) represents the instantaneous variation of N . The second term is the vertical advection of droplets: w is the vertical wind (component of the 3-D-wind W along the vertical axis u z ); v is the terminal velocity of droplets under gravity. In this study, the parameter v is assumed to depend only on the diameter D, according to the Atlas et al. (1973) with D in m and v in m s −1 . This relation was fitted to be in good agreement with Gunn and Kinzer (1949) measurements (which are often used as a reference) in the diameter range [0.6-5.8 mm] (Atlas et al., 1973). Equation (1) can be considered as a simplified version of the model proposed for instance in Hu and Srivastava (1995). We will now present this complete model and discuss its different terms.
2.1.2 Complete propagation model (Hu and Srivastava, 1995) and discussion of its different terms In the (Hu and Srivastava, 1995) study, N depends on time, three coordinates of space and diameter (N = N(t, x, y, z, D)) and obey the PDE: We recognize the two first terms of Eq. (1) with the difference that we have a 3-D wind W advection here. The third term of Eq. (3) represents the evaporation. τ is the net rate of change of droplets' masses (unit kg s −1 ) and dD/dm is the derivative of the diameter-mass relation for spheric drops (unit m kg −1 ). The last term I represents the collision effects (drop mass changes not due to evaporation).

3-D wind W
The vertical wind is considered to add an offset to the terminal velocity of drops. Additionally, the knowledge of the real droplet vertical speed is critical to have information about the DSD from the return power spectra of a vertically pointing Doppler radar. This question has been widely investigated in the literature (Lhermitte, 1988;Giangrande et al., 2012;Tridon et al., 2013) and it shows that we have to take the vertical wind into account to properly deal with Doppler radar data. The horizontal wind is the main force which makes the droplets move in the horizontal plane. However, with only one vertically pointing radar and co-localized ground point measurements of the DSD, it is not possible to record the horizontal variability of the DSD. Consequently, as mentioned above, we limit the study to a (z, t) plane. This simplification is consistent if the horizontal air motion is weak or if the DSD is homogeneous in the horizontal plane. (Of course the DSD can be inhomogeneous with height.) To limit the errors caused by this simplification, we will only treat quiet, stratiform light rain events (see the case study in Sect. 3.3) and we will study the ability of the algorithm to reproduce the observations.

Evaporation (third term of Eq. 3)
The evaporation is mathematically modeled as a term of advection along the mass of droplets coordinate. The advection speed (parameter τ of Eq. 3) depends on the drop diameter, and on several meteorological variables, including pressure, temperature and mainly humidity. Seifert (2008) proposes a parameterization of this term.
According to this parameterization and for long cold stratiform rain events similar to the one investigated in this study (see Sect. 3.3), we verified (not shown) that evaporation can be neglected.

Collisions (fourth term of Eq. 3)
Collisions between drops can lead to coalescence (two drops producing one drop) or breakup (two drops producing at least two drops). Coalescence/breakup phenomena are generally assumed to be an important factor governing the DSD temporal evolution (Mcfarquhar, 2010;Prat et al., 2012), even if some studies reconsider this assumption (Villermaux and Bossa, 2009;Villermaux and Eloi, 2011). The phenomena have been widely investigated. Numerous parameterizations have been proposed. Among them, some characterize the ability of two drops to coalesce, depending on the energy involved in the collision (Low and List, 1982b;Brazier-Smith et al., 1972. Others are more interested in the distribution of resulting drops in the case of breakup. The resulting distributions can be based on laboratory experiments, Low and List, 1982a, b, or on numerical fluid dynamics models, Schlottke et al., 2010; Since our main objective is to retrieve DSD profiles using data assimilation, we restrict ourselves to a simple framework, and we do not take into account the collisions between drops. Moreover, among the wide literature mentioned above, many authors show that the coalescence/breakup processes are less critical for low rain rates (List et al., 1987;Prat and Barros, 2007a;Barthes and Mallet, 2013). Thus, this paper will only focus on a light rain event in order to limit the errors caused by this simplification. Additionally, as for the horizontal wind exclusion, we will study the ability of the algorithm to reproduce the observations despite this simplification.

Numerical model (discretization)
This simplified PDE used in this study (Eq. 1) is discretized on an (time/height/diameter) Arakawa C-grid (Arakawa and Lamb, 1977), meaning that N is evaluated at the grid-box centers, while w is evaluated at the grid-box faces. For this discretization of the PDE, we use the finite difference scheme of Smolarkiewicz (1983), developed for advective problems in atmospheric modeling requiring a large number of variables to be dealt with. For this scheme, the Courant-Friedrichs-Lewy stability condition is with t and z the time and height steps, respectively. In order to be consistent with radar observations (see Sect. 3.1), we choose z = 100 m. Then, according to the usual range of drop terminal velocities (v ∈ [0; 10] m s −1 ) and of stratiform vertical wind speeds (w ∈ ±2 m s −1 ), we choose t = 5 s. The number of time and height discretization steps are denoted NT and NZ, respectively. In order to have a good coverage of the range of diameters, they are discretized from D min = 0.2 mm to D max = 7.5 mm, with a diameter step D = 0.1 mm. The resulting number of diameter bins (ND) is thus equal to 73. We note N = ( N , the discretized DSD field. N n i indicates a specific discretized DSD at grid point (n, i), namely, at time t n = t 0 +n t and height z i = z 0 −i z. Similarly, the notations w and v stand respectively for the discretized wind field and terminal velocity vector.
The unknowns of such a model are as follows.
1. The initial (t 0 ) DSD field (initial condition), , is unknown. However, if we suppose that the model starts running before the beginning of the rain event, we can suppose that this initial field is 0 everywhere. This is the assumption made hereafter.
3. The vertical wind time-height field, w (which is, in our case, a parameter of the model), represents NT × NZ unknowns.
In order to reduce the dimension of the problem, and thus to substantially reduce the number of degrees of freedom, we have to add some a priori information. Since the density of raindrops in the atmosphere is well approximated by a gamma distribution (Mcfarquhar, 2010), we parameterize the top DSD ( N with α the total volumetric number of droplets (unit m −3 ), and f the gamma probability density function (unit m −1 ): with k the shape coefficient and θ the scale coefficient. Thanks to this parameterization, for each time step n, we limit the top DSD ( N n 0 ) n∈ [[1:NT]] determination to the choice of 3 parameters and so divide the resulting number of unknowns by almost 25 (from 73 to 3). We note x, the corresponding three parameters' field (dimension: 3 × NT): . Concerning the implementation of the algorithm, at each time step n, x n is directly converted into a top DSD N n 0 using Eqs. (5) and (6); therefore the propagation model can use it as input data.
Finally, the problem is reduced to the determination of the parameter field ( x) along with the vertical wind field ( w), resulting in a total number of unknowns of (NZ + 3) × NT.

General principles
In the previous section, we have seen that the model propagates the top DSD ( N n 0 ) n∈[[1:NT]] through space and time according to the vertical wind w. Besides the numerical model, we have a full set of observations at various heights and times (from radar and disdrometer; see Sects. 3.1 and 3.2).
The retrieval of the top DSD and of the vertical wind is thus carried out through the use of variational data assimilation which basically consists in minimizing the distance between the model and the observations (see Fig. 1 which details our contextualized assimilation process).
To perform a 4-D-VAR assimilation, we still have to build a 4-D-VAR cost function, which is crucial since it is used to evaluate the gap between the observations and the top DSD ( N n 0 ) n∈ [[1:NT]] propagated by the model to different times and heights. The numerical minimization of the cost function allows a set of unknowns to be achieved (namely, x and w) which minimize the cost function.
In the data assimilation context, the cost function usually takes the form: with N the discretized DSD as previously defined. J o is the observation term of the cost function, presented later. J b is the background term (also called first guess) of the cost function. J b keeps the solution in the vicinity of a given a priori state. In our case, we suppose that there is no background available; therefore J b will be dropped from the cost function. Moreover, we add regularization and penalization terms (see Sect. 2.2.4). The cost function finally used takes the form: where the cost function is the sum of a regularization term (J r ) along with two penalization terms (J x ), dedicated to ( x) the DSD parameter field, and J w , dedicated to ( w) the wind field. These three terms will be discussed and explained in Sect. 2.2.4. In the end, the cost function J ( x, w) is a function of both the top DSD parameter field ( x) and the wind field ( w) that allows the overall distance between a propagated set of parameters ( x, w) and the observations to be computed. By minimizing the cost function, the algorithm will estimate an optimal set of ( x, w), minimizing this distance. Once the assimilation process is completed, for each time step, we will have a top DSD parameter and wind values for all heights.

Observation term of the cost function J o
J o , the observation term of the cost function, takes the form: For a given grid point (n, i), y n i is the vector composed of the observations from radar and disdrometer. H i,n is the operator projecting the DSD on the observation space. Note that, here, the covariance matrix of the observation error is the identity.
At time (t n ), using the vertical wind field ( w), the model has propagated the top DSD N n 0 from times t n < t n . We note M n i the numerical model propagating all the necessary top DSD N n 0 to compute N n i , the DSD at time (t n ) and height (z i ). Then, we have with x the three DSD parameters field as defined above (see Sect. 2.1.3). Consequently, the cost function (J o ) can be expressed as a function of the unknowns x and w only:

Observation operator H
The observation operator, H, can be considered as an aggregation of 2 sub-operators H dis and H dop , modeling the disdrometer and radar observations, respectively. For disdrometer measurements, we build an operator H dis able to convert the propagation model ground into disdrometer-like measurements (i.e., drop histograms in our case). A disdrometer records a flux of droplets crossing a given area (A, in m 2 ). Then, at time t n and for drops in the diameter class p, this operator has the form: For radar measurements, we build an operator H dop converting model DSD into Doppler radar-like measurements.
Atmos. Meas. Tech., 9, 3145-3163, 2016 www.atmos-meas-tech.net/9/3145/2016/ Radars operating in a vertically pointing mode provide the distribution of the radar return power as a function of the Doppler velocity, the so-called Doppler spectrum (Giangrande et al., 2010). The retrieval of the DSD directly from these Doppler spectra is an inverse problem (Kollias et al., 2002) mixing various phenomena. Here we build a "direct" operator which acts in reverse order (from DSD to Doppler spectra). It takes vertical air motion and rain attenuation into account. At a given time and location, this operator can be computed according to the three following steps.
1. Calculate the theoretical non attenuated and non-noisy spectrum F (v) (s m −2 ) from a DSD N (D): where, given a drop diameter D, σ (D) is the backscattering cross section (unit m 2 ), calculated for the radar frequency according to the Mie scattering theory (Ulaby et al., 1982). dD/dv is the derivative of the diametervelocity relation.
2. Then derive the final attenuated spectrum F a (v), encompassing both the effects of the attenuation and of the vertical air motion w.
The exponential part of Eq. (14) adds the effect of rain attenuation (see for instance Peters et al., 2005), with K(z) the attenuation coefficient at height z (unit m −1 ), calculated through K(z) = D σ ext (D)N (D, z)dD, with σ ext (D) the extinction cross section (Mie theory).
3. Finally, for minimization purposes, F a is converted into a logarithmic scale, resulting in a logarithmic bounded attenuated spectrum Z(v): We recall that the operator run on the computer is necessarily a discretized version of what is presented here.

Regularization and penalization terms
In this section, we present a brief description of the penalization term J r as well as the two regularization terms J x and J w of the cost function (see Eq. 8).
1. J r : we have seen in Sect. 2.1.3 that the time step of the model is set to 5 s. It means that x and w have to be retrieved at this resolution. However, the minimal resolution of the instruments is generally greater than this value (10 s in our case; see Sects. 3.1 and 3.2). Then, using J r avoids the assimilation algorithm producing highfrequency fluctuations on both of the unknowns x and w by giving lower cost to smooth fields.
2. J x allows us to proceed to the minimization under constraints for x, the three DSD parameters penalization term. Because the instruments used are not able to record the complete drops diameter range, it is used to avoid the assimilation algorithm producing unphysical DSD. To this end, we add this penalization term which avoids the components of the x DSD parameters field to go out of predefined ranges: [0; 8000] m −3 for α; [0; 3] for k; [0; 10 −3 ] m for θ .
3. J w : the vertical wind field w modifies the advection velocities and shifts the Doppler spectra. Then, the vertical wind is mainly controlled (in the assimilation algorithm) by the Doppler spectra. If no spectra are available for a given vertical layer, the assimilation algorithm could produce unphysical temporal fluctuations of the vertical wind which would not imply significant extra cost without the term J w . Since there may be no Doppler spectra available at given heights (see Sect. 3.1), we use J w , the wind penalization term, to force the vertical wind to stay close to 0 by adding extra cost to strong vertical winds.

Minimization of the cost function and parameter estimation
The numerical minimization of the cost function J ( x, w) goes through the calculation of its gradient. Note that all the regularization terms are easily differentiable. The gradient of the observation term (J o ) of the cost function is with t M n i the adjoint of the linearized version of the operator M n i . We use the YAO software (Thiria et al., 2006;Nardi et al., 2009), developed by the LOCEAN (Laboratoire d'Océanographie et du Climat). It provides a simple method for deriving the adjoint of the model. For the numerical minimization itself, it is coupled with M1QN3 (Gilbert and Lemaréchal, 2006), a quasi-Newton algorithm to solve unconstrained optimization problems. Once the propagation model is implemented, it allows the use of observations to compute realistic estimations for the unknowns.
We finally note that we have to initialize the unknowns of the problem ( x, w) at the very beginning of the assimilation process. We found that this initialization is not critical. Consequently, we will initialize the wind field ( w) at 0 and the 3 gamma parameters time series ( x) at constant very low values, namely α = 1, k = 0.8, θ = 2.10 −4 .
In this section, we firstly give the characteristics of the two instruments used in this study, namely a Doppler 24 GHz radar, and a disdrometer. Then, we present the characteristics of the case study examined.

Micro-rain radar (MRR)
The first instrument used in this study is a micro-rain radar (MRR, Peters et al., 2005), developed by METEK GmbH, Germany, and belonging to Météo France. The MRR is a 24 GHz frequency modulated, continuous wave (FM-CW), vertically pointing Doppler radar, with a small transmit power (50 mW) and a beamwidth of 2 • . The raw return power spectrum is processed by the MRR (see METEK (1998)), which supplies the spectral reflectivities. In our assimilation process, only the spectral reflectivities will be used as observations (noted F MRR ) .
The MRR software also supplies DSD estimations (obtained through an inversion algorithm of the Doppler spectra, including attenuation correction) and estimations of different moments of the DSD (liquid water content, rain rate, reflectivity factor). Those quantities, mentioned as "MRR products", will be used as a comparison framework for the DSD N retrieved by our algorithm.
The radar range gate size is 100 m. The MRR provides data up to a 3100 m height. The temporal resolution is 10 s. This fine temporal resolution will allow us to integrate the observations over various time periods. This integration time will be referred to later as an "observation window". The Doppler velocities used a range from 0.56 up to 9.54 m s −1 , with a step of 0.19 m s −1 . According to the manufacturer recommendations (METEK, 1998), the MRR records for the two first radar gates (< 200 m) will not be used in the assimilation process.

Dual-beam spectropluviometer (DBS)
The DBS is an optical disdrometer developed by the LAT-MOS, whose main advantage is its ability to cover a very large range of raindrop diameters . Its collecting area is 0.01 m 2 . For each recorded drop, the DBS supplies a triplet (time of arrival/spherical equivalent diameter/fall velocity), allowing the measurement of drop flux time series, noted N DBS , used as observations in the assimilation algorithm. We apply a minimum diameter threshold of 0.4 mm. 0.4 mm corresponds to a value over which we are sure that the instrument avoids false detections  and can be compared to other instruments (Krajewski et al., 2006).We will only conserve the drops above this diameter.

Case study data set
The two instruments described above were deployed in Ardèche (southwest of France), during autumn 2013, in the context of the HyMeX campaign (HYdrological cycle in the Mediterranean EXperiment, see for instance Bousquet et al., 2014). Figure 2 presents a topographic map of the area. The two instruments (DBS and MRR) were co-localized at Le Pradel (44.6 • N, 4.5 • E), in a mountainous area called the Cévennes-Vivarais, subject to so-called Cévenol flash flood events (Molinié et al., 2012). However, for the purpose described in Sect. 2.1.2, we will not work on extreme rain events, but rather on a stratiform event of light rain which occurred on 12 October 2013.
This event occurred at Le Pradel in the afternoon, from 16:30 to 21:30 UTC. It is a long rain event characterized by low but sustained rain rates at ground level (the disdrometer was used to record 2 min rain rates always between 1 and 6 mm h −1 ; see Fig. 3, red line). The event is also quite homogeneous. All the Météo France automatic weather stations located around Le Pradel (see Fig. 2) record consistent cumulative rainfall: 11.5 mm at Le Pradel, 12 mm at Aubenas, and 12.9 mm at Berzème. From the MRR (Fig. 4, left) and DBS (Fig. 3, red) measurements, we can distinguish two distinct phases in the event. Until 18:15, there is very light rain, with only small drops (rain rates around 2 mm h −1 and very low reflectivities and characteristic velocities). After 18:15, there are short periods with higher rain rates (4-5 mm h −1 ) and, mainly, higher reflectivities and falling velocities (up to 35 dBZ and 8 m s −1 ).
The temperature was around 8.5 • C at ground level (Le Pradel station, 278 m height). Assuming a gradient of −1 • C/150 m, we found a freezing level at 1550 m, which is consistent with the MRR records (on Fig. 4, left, we see the so-called radar bright band, indicative of melting ice crystals, at about 1500 m height). As a consequence, the top of the retrieval domain in our assimilation algorithm will be set to 1400 m in order to avoid bright band problems (the propagation model is not able to treat the melting layer). Note Atmos. Meas. Tech., 9, 3145-3163, 2016 www.atmos-meas-tech.net/9/3145/2016/ For this case study, the number of time and height discretization steps (NT and NZ; see Sect. 2.1.3) is set to 4900 and 15, respectively. We remind the reader that, for numerical purposes, the time and height steps were chosen to be t = 10 s and z = 100 m.

Simulated data (twin experiment)
We have seen in the last section that the case study chosen to apply our retrieval algorithm is consistent with the propagation model underlying hypothesis. As a reminder, it assumes DSD to be driven by gravity and vertical wind alone. Thus the effects of evaporation, collisions and horizontal wind are discarded. However, other phenomena remain unaccounted for in the model. Since this could impair the ability of the model to satisfactorily manage our real case study, we have to assess the impact of two of these phenomena (instrumental noise and turbulence) on the assimilation retrievals. This is done, in this section, by performing a twin experiment, whose principle is described below.

Description of the twin experiment process
Twin experiments consist in applying the assimilation process to simulated data with properties analogous to real data. Figure 5 schematically summarizes all the successive steps of our twin experiment (i.e., to simulate realistic observations and to run the assimilation). Firstly, we simulate realistic series of top DSD parameters and vertical winds (step a on Fig. 5). These parameters are what we would like to retrieve with the algorithm. In Sect. 4.2, we will explain how these fields are simulated. Then, these series are propagated through space and time using the propagation model (step b) and "true" (namely, by supposing the model to be perfect) DBS and MRR data are produced through the use of the observation operators (step c). All the fields obtained through this process (DSD and vertical wind, as well as observations, for instance Doppler spectra) are mentioned below as true.
Then, from these fields, we simulate realistic MRR and DBS observations by adding model and instrumental errors (instrumental noise and turbulence, see Sect. 4.3), which cannot be directly produced by the propagation model (step d) since they are not present in the model despite their effects on the measurements. These data are mentioned below as "observed". The process used to simulate these observed fields is detailed in Sect. 4.3.
Finally, the assimilation algorithm estimates an optimal set of unknowns explaining these observations (step e). Assimilated x and w fields are then propagated by the model (step f) and the observation operators (step g). All these retrieved fields are mentioned below as "assimilated". They will be compared to true and observed fields. The results of this data assimilation twin experiment are detailed in Sect. 4.4.
This twin experiment demonstrates the ability of the assimilation algorithm to produce coherent DSD fields when the assumptions of our propagation model are satisfied. It also shows how it indirectly handles the turbulence, underwww.atmos-meas-tech.net/9/3145/2016/ Atmos. Meas. Tech., 9, 3145-3163, 2016 lining its efficiency despite the presence of non-modeled phenomena. Note that in this section, all the observations (MRR spectra F MRR and DBS DSD N DBS ) are integrated over a 2 min observation window (if not, it will be explicitly specified).

Simulation of top DSD parameters ( x) and vertical wind field ( w)
First, we have to simulate realistic time series of top DSD parameters ( x). To do this, we use the rain rates recorded on the ground by the DBS during the case study presented above (Sect. 3.3). From each of the successive disdrometer rain rates, we calculate the two parameters of an exponential DSD according to Marshall and Palmer (1948). This time series of Marshall-Palmer DSD (special case of gamma DSD) is our simulated top DSD ( N n 0 ) n∈[[1:NT]] , which is used as input in the propagation model.
We also need a time-height field of vertical winds w, (resolution 5 s, 100 m) to run the propagation model and simulate observations. We have seen that our algorithm is, up to now, limited to stratiform quiet rain events. For this kind of event, vertical winds range in ±2 m s −1 , with temporal correlations around 2-5 min and spatial correlations around 500-1000 m (Giangrande et al., 2010). The process used to simulate such a field is as follows. (1) We create a 5 s/100 m field of independent, normally distributed, random variables. (2) We average this field over a mobile (2 min/500 m) window. (3) We scale the field to get winds within ±2 m s −1 . In the end, we obtain a wind field with triangular correlations corresponding to the expected values. Figure 7 (top, on the left) shows a part of the vertical wind field used in the twin experiment (i.e., on simulated data). Note that the winds are positive downward.

Simulation of realistic observations
We detail the process used to simulate realistic observations in this section (addition of turbulence effects and instrumen-Atmos. Meas. Tech., 9, 3145-3163, 2016 www.atmos-meas-tech.net/9/3145/2016/ tal noise). The purpose of this step is to build observations consistent with what is expected from real measurements.

Radar Doppler spectra
We want to simulate Doppler spectra as they would be recorded by an MRR in realistic conditions. For this purpose, we run the propagation model and save the attenuated reflectivities calculated by the algorithm (observation operator) from the DSD (variable F a (v), see Eq. 14). Then, for each spectrum, we apply the following process.
(1) The atmospheric turbulence impact is added. Turbulence can be seen as a fine-scale modulation of the vertical wind. For a time and location, the effective vertical air motion is the sum of the mean vertical wind (defined above and taken into account in the propagation model) and of a random variable, modeling the turbulence (not taken into account). Mathematically, the impact of turbulence on the Doppler spectra can be modeled as the convolution of the Doppler spectra with a Gaussian function (Gossard, 1994;Williams, 2002;Tridon et al., 2013). So, the turbulence is added to our spectra by applying a discretized version of the convolution operator: with * the convolution operator and σ t the width of the turbulence spectrum. We use σ t = 0.7 m s −1 , considered a realistic but large value (Peters et al., 2005;Tridon et al., 2013) implying the use of a large spectral broadening ).
(2) Gaussian white noise (intensity 1 dB) is added to F t . On Fig. 6, on the right, the dotted lines show observed Doppler spectra at different heights.

Disdrometer
We simulate the records of a ground disdrometer. For this purpose, we run the propagation model and save the discretized DSD in the grid box just above the ground. For a given time and diameter class p, we note n p , the corresponding number of drops in a given box. Then, for each of these n p drops, we choose a diameter (following a uniform law in the interval [p; p + 1] D + D 0 ) and a velocity. This velocity is the theoretical Atlas velocity for the diameter modulated by Gaussian white noise (as presented previously), representing the atmospheric turbulence. The height of the drop is finally uniformly drawn in the interval [0; 100] m (height range of the grid box). At the end of the process, we check which drops, among the n p , will touch the ground during the time interval t (and so will be seen by the disdrometer). The drops touching the ground are recorded as observed by the disdrometer. This process is a way to mimic the natural variability of rain and to run the observation operator defined by Eq. (12) in a stochastic way.

Results -robustness of the algorithm
In this section, we firstly assess the ability of the assimilation process to reproduce the observations. Then, we compare the true and assimilated DSD N and wind fields w. Finally, we investigate how the model response is impacted by the size of the observations window. First of all, we have to introduce the numerical indicators needed to perform this investigation.

Error indicators
In this section, we define the indicators used to quantitatively evaluate the performance of the algorithm. Choosing error in-dicators is rather sensitive to the range of values under study. For fields with magnitude far from 0 and ranging in small intervals, we use relative indicators without any risk of giving too much importance to very low (close to 0) data. We define the two following indicators: the mean absolute percentage error (MAPE), and the relative bias (rbias), defined by where A and B are the respective reference and evaluated fields rearranged in vector forms. P stands for the dimension of the vectors. Fields with very low or negative values render the two previous indicators inappropriate. For these fields, we use the mean absolute error (MAE), defined by

Observations reproduction
In the following analysis, we use the denominations true, observed and assimilated (see Sect. 4.1 and Fig. 5). Table 1 gives the MAPE and the relative biases between true, observed and assimilated characteristic velocities and spectral widths. We recall (Sect. 4.3.1) that the atmospheric turbulence, which broadens the Doppler spectra, is included in the observations but not in the true data. Consequently, for the observed spectra, we obtain spectral widths larger (+12 %) than the true ones.
The assimilated spectra are also broader than the true ones, but to a lesser extent (spectral width 9.5 % larger). They are slightly less broad than the observed spectra (−2.3 %). We also note that the MAPE and rbias have almost the same absolute values (Table 1, third and fourth rows). This means that the sign of the differences between the two compared fields is the same for almost all the grid points. This implies that all spectra are broadened by turbulence independently of time and height.
On the contrary, the characteristic velocity is not highly affected by the addition of turbulence impact on observed fields (bias −0.36 %, MAPE 0.40 %). Consequently, the characteristic velocities are correctly reproduced by the algorithm (biases and MAPE inferior to 2 %). Obtaining the result that turbulent broadening only impacts the spectral width and not the reflectivity-weighted mean downward velocity in the simulations confirms the robustness of the assimilation procedure and expected results from prior work (Williams, 2002;Williams and Gage, 2009).
To reproduce the turbulence broadening, the algorithm can act on the two unknown fields, x and w. Firstly, it can affect the DSD shape, by adding small and large droplets and removing drops of intermediate diameters; but because the DSD is directly controlled on the ground by the disdrometer (observation N DBS ), it forces the algorithm to reproduce the ground DSD ( N n NZ ) n∈ [[1:NT]] well. Between true and assimilated disdrometer observations, we find an MAE (see Table 1, last row) of 27 mm −1 , while it is 42 mm −1 between truth and observations (the mean number of drops being 689 mm −1 ). Secondly, the algorithm can use another way, which consists in using the fast modulations of the vertical wind w at fine resolution. We have seen that all the observations are integrated over a 2 min observation window (Sect. 3.1), while the wind is retrieved at the 5 s model resolution. Consequently, if the 5 s winds oscillate inside the 2 min observation window, it will result in shifting the 5 s corresponding spectra alternatively towards small and large Doppler velocities. Once averaged over the 2 min time period, these shifted spectra will result in a final 2 min spectrum with a spectral width larger than all the 5 s spectral widths, whatever the original DSDs. This is the way used by the algorithm to mimic the turbulence. This phenomenon can be seen in the third column of Table 2, which shows the MAE between true and assimilated wind fields for an observation window of 2 min and for various field resolutions (spreading from 10 s to 8 min). We can see that the MAE of the 10 s resolution wind fields is large (0.46 m s −1 ) compared to the true absolute mean value (0.34 m s −1 ; see last row of Table 2), but it sharply decreases for resolutions greater than 2 min. This means that there are fluctuations at the 5 s model resolution which do not correspond to the physics. These are numerical artifacts produced by the model which uses any degree of freedom available to reproduce the broadening observed in the measurements.
An example of Doppler spectra at a 2 min resolution is given in the right part of Fig. 6. In this plot, examples of vertical profiles of Doppler spectra are shown. On the observed and assimilated spectra, we can clearly see the turbulence broadening, absent in true ones. Since this broadening on assimilated spectra is due to a turbulence reproduction (using the wind resolution) and not to a modification of the DSD, the good behavior of the assimilation algorithm for realistic stratiform rain is thus assessed.

Wind and DSD reproduction
Even if we confirmed the ability of the assimilation algorithm to handle the turbulence, we still have to assess some other important rain properties.
As in the previous section, we compare the MAE between true and assimilated fields for an observation window of 2 min and for various resolutions in Table 2. The quantities Z e and D m (moments of the DSD, see caption of Table 2) are very well reproduced by the assimilation algorithm, with very small MAE relative to the true absolute mean values Atmos. Meas. Tech., 9,[3145][3146][3147][3148][3149][3150][3151][3152][3153][3154][3155][3156][3157][3158][3159][3160][3161][3162][3163]2016 www.atmos-meas-tech.net/9/3145/2016/  Table 2. For the twin experiment: MAE (Eq. 20) between true and assimilated fields for different parameters: vertical wind; rain rate, RR; reflectivity factor, Z e ; mean volume diameter, D m ; liquid water content, LWC; number of drops, N 0 . The indicators are presented for two observation windows (10 s, in italics, and 2 min) and for different temporal field resolutions (from 10 s to 8 min). The absolute bias is also presented, as well as the absolute values of the means of the true fields. In each case, the best result (lowest MAE) has been made bold if it happened in the case of the 2 min observations window.

Field
Obs (for instance, at 2 min field resolution, MAE of 0.80 dB and 24.7 µm compared to true absolute mean values of 27.6 dBZ and 1.12 mm, respectively). The results are slightly less satisfactory for liquid water content (LWC) (10.0 vs. 142 µg m −3 ) and for N 0 (241 vs. 1070). This is due to the fact that the number of drops is mainly driven by the number of very small drops (< 0.4 mm), which are not seen by the disdrometer. Moreover, given their very low backscattering cross section, they have almost no weight in the cost function. They are consequently not effectively controlled. However, we can note that these drops represent a very low mass of water and that their non-control is not critical. Moreover, since the top DSD is parameterized, we cannot obtain unreasonable values.
We have seen that the 2 min MAE for the wind field is quite satisfactory, meaning that the global trends of the wind field have been reproduced by the algorithm; but we would also like to know whether the wind patterns are well reproduced. Figure 7 shows extracts of the true and assimilated wind fields, as well as the mean and standard deviations of these fields according to the height (radar gate) and the temporal autocorrelation functions. The top images confirm a satisfactory reproduction of the wind fields. We also see that the correlation time, the absolute mean values and the standard deviations are well reproduced by the assimilation, even if the standard deviation results are more dubious for the two lowest radar gates, where no Doppler spectra are available (as seen in Sect. 3.1).
The rain rates are also satisfactory for resolutions over 2 min (Table 2, column 4), with an MAE of 0.14 mm h −1 with a mean value of 2.18 mm h −1 .

Impact of the observations window
In this last subsection, we briefly assess the influence of the integration time of the observations (observation window, see Sect. 3.1). Table 2 presents the MAE between different fields, depending on the observations window, 10 s or 2 min. We have seen that a 2 min observations window allows the algorithm to use the 5 s model resolution to reproduce the turbulence spectral broadening by adding artificial synthetic fluctuations. If the observation window is reduced to 10 s, no room is left for fluctuation. Consequently, the MAE for the wind is smaller with 10 s observations window for fine resolutions (for instance, 0.21 vs. 0.46 m s −1 at 10 s).
Because the model fits, at best, the 10 s resolution, it is natural that the algorithm performs less accurately at the 2 min resolution and hence beyond.
The results are globally the same for the other parameters of Table 2. Because the algorithm cannot use the wind w to adjust the assimilated spectra to the observed broadened ones, it has no choice but to fit the observations through a change of the top DSD (excess of small and large drops) in order to produce the best compromise at the cost level. This effect is visible in the N 0 column of Table 2: the results are worse for a 10 s observation window whatever the resolution (excess of small droplets). We can also note that the spectra are much noisier at 10 s, implying a tougher minimization.
Finally, we conclude that an observation window of 2 min is a better choice than 10 s. It allows the algorithm to handle the turbulence and limit the instrumental noise influence. Thus the assimilation in the next section will be run with a 2 min observation window.

Results for real case study
In this section, we apply the assimilation algorithm to the real data case study described in Sect. 3.3. Again, we show that the assimilation algorithm is efficient (able to fuse different observations) for this case study. We then show that it produces spatiotemporally coherent wind and DSD fields.

Observation reproduction
We use the indicators introduced in Sect. 4.4.1 again (see Table 3). Firstly, we note that the disdrometer observations N DBS are well reproduced. The MAE between observed and assimilated data is 102 mm −1 , while the mean observed number of drops is 1080 mm −1 . These results are very similar to the results of the previous simulated study, meaning that the algorithm successfully reproduces the disdrometer observations. Figure 3 shows the ground rain rates as calculated from the disdrometer observations (red) and as reproduced by the assimilation algorithm (blue). It confirms that the algorithm is efficient for ground data. We will later come back to the remaining differences between the two rain rate time series. Table 3 also presents results concerning the MRR observations. Even if the MAPEs for the characteristic velocity W c and spectral width SW are larger than the ones on simulated data, they remain satisfactory (mean error of 4.4 and 14 %, respectively).
The three columns on the right of Table 3 present these results at different heights and help us to analyze the vertical structure of MRR observations and assimilated data. When considering reflectivity η (zeroth order moment, second row), we see that there is globally no bias between observed and assimilated reflectivities (−0.017 dB). However, the assimilation process overestimates the reflectivities at the top (+0.60 dB) and underestimates the reflectivities close to the ground (−0.40 dB), while there is a smaller bias at middle heights (−0.12 dB). Since the model propagates the drops from the top to the bottom, without changing their diameter, the only process in our propagation model that can produce a decrease of the mean reflectivity with height at the scale of the event is attenuation. This means that the MRR observations display more reflectivity decreasing with height than what can be explained by the model attenuation. This feature is confirmed by the left part of Fig. 4. We clearly see that the rain rate (top) and reflectivity factor (middle) display a decrease with height for the MRR products. This phenomenon cannot be reproduced by our algorithm (right part of Fig. 4) which more or less (depending on the wind field, see Sect. 5.3 below) conserves these quantities with height. We also see that the algorithm reaches a compromise between the different heights; therefore the assimilated reflectivities (and hence the rain rates and other DSD parameters displayed in Fig. 4) are overall (compared to MRR products) slightly overestimated on the top and underestimated on the ground. This would also explain why, in Fig. 3, assimilated rain rate peaks are generally underestimated compared to disdrometer data. Figure 8 shows the bottom DSD (continuous lines) for the MRR products (200-300 m, red), for the DBS measurements (on the ground, in black), as well as the assimilated ones (on the ground, in blue). For the intermediate diameters ([0.7; 2.5] mm, corresponding to the range which mainly impacts the rain rate and the reflectivity factor at low rain rates), the assimilation underestimates the number of drops. This could be explained by the MRR observed spectra at the top, which induce a lower number of drops (Fig. 8, red dotted marked line). The origin of this variation of reflectivity with height for the MRR observations is not clear. It could be a calibration problem (heightdependent transfer function; see METEK, 1998), as well as the effects of the horizontal wind or of some microphysical Atmos. Meas. Tech., 9, 3145-3163, 2016 www.atmos-meas-tech.net/9/3145/2016/ processes (coalescence/breakup) assumed negligible in our propagation model. Anyway, the assimilation algorithm appears able to merge the different data available to produce a solution, making a compromise between the observations available (low biases and MAE/MAPE for the three first moments of the spectra and for disdrometer data; see Table 3). It produces results embedding the spatiotemporal coherence brought by the propagation model. We can also note that the very small drops (< 0.4 mm) are loosely constrained and that the corresponding results are then dubious (Fig. 8). The disdrometer does not cover these diameters and the MRR does not receive a lot of energy for the corresponding Doppler velocities, so that the MRR observed spectra seem less coherent. Consequently, the assimilation algorithm produces a very large number of small drops (more than an exponential law would) without any physical justification associated.
For large drops (> 2.5 mm), we note (Fig. 8) that the disdrometer provides observed ground DSD N DBS which are not reproduced by the assimilation, probably because they are not consistent with the MRR observed spectra. However, large drops are very rare for this light rain event; therefore the presence or absence of a particular drop can deeply affect the tail of the disdrometer DSD. Moreover, the terminal velocity of these drops is > 7.3 m s −1 ; therefore disdrometer measurements with 2 min observation window correspond to very large probed volumes (corresponding to a height of around 850 m, crossed in 2 min for a 7.3 m s −1 drop) which absolutely do not match the MRR probed volumes (100 m height). Because large drops are rare for this event, we nevertheless chose to work with a 2 min observation window. Naturally, for convective events, this choice could be questioned, at least for the disdrometer (we have seen in Sect. 4.4.4 that a decrease of the observation window width for the Doppler spectra decreases the ability of this algorithm to handle the turbulence). We remind that there are no Doppler spectra available for the two lowest gates (0-200 m); therefore the wind retrievals for these two heights are not constrained by observations in the algorithm, since the wind also acts on the vertical advection of drops and therefore also on the disdrometer records.

Wind
The structure of the assimilated wind field is consistent with what we could expect for such an event. The correlation characteristic time and height are around 1 min and 300 m for the 10 s wind field (which can be used by the algorithm to handle the turbulence as explained in Sect. 4.4.4) and around 3 min and 500 m for the 2 min field. These results, as well as the wind range (±1.5 m s −1 ), are consistent with the results of Kollias et al. (2002) or Giangrande et al. (2010), for instance. Even if we cannot validate the retrieved field using independent measurements in this study, we could nevertheless notice that our method produces the same wind field structures whatever the algorithm initialization.
We also have an overall mean wind of +0.004 m s −1 , slight and not notable downward air motion. We get slight downward air motion in the upper layers of the atmosphere (mean wind around +0.15 m s −1 above 700 m) and slight upward air motion closer to the ground (around −0.2 m s −1 between 200 and 500 m). The 2 min standard deviation is quite constant with height and around 0.25 m s −1 .

DSD
We have seen (Sect. 5.1) that, on average, the MRR products (rain rate, reflectivity factor) are consistent with our assimilated data. There is nevertheless a vertical decreasing tendency on MRR products, which cannot be reproduced by the algorithm, given that the propagation model produces verti-cally coherent fields (Fig. 4). However, in this section, we focus on the structure of the assimilated DSD fields, independently of the MRR products.
We have seen that the algorithm conserves the reflectivity factor during a drop fall well. It results in diagonal structures on the time-height plane (see Fig. 4, middle, right). The slope of a diagonal structure reflects the fall velocity. We can see that the reflectivity field is consistent with the characteristic velocity field (bottom, right): the lower the velocity, the lower the slope of the reflectivity diagonal structure. Now we will also assess whether these structures are retrieved for some important moments of the DSD. Figure 9 shows the mean volume diameter D m , and the number of drops N 0 , resulting from the assimilated DSD N . We note that the retrieved values for D m (0.5-1.5 mm) are reasonable for light stratiform rain (see Fig. A4 in Peters et al., 2005). D m is also quite well conserved with height, but the patterns seem noisier than for the reflectivity, as it also does for the rain rate field (Fig. 4). This is a consequence of the vertical wind, which modifies the vertical advection and so the DSD. The vertical wind has a greater impact, relatively, on small drops than on large ones (given their smaller terminal velocity); so the low-order moments, mainly driven by the small drops, are more affected than the high-order moments. This is particularly clear for N 0 . We have seen that the algorithm, underconstrained for these diameters, produces a very large number of very small drops (< 0.4 mm, i.e., with terminal velocities < 1.5 m s −1 ) (see Fig. 8). Consequently, (1) the vertical consistency of N 0 is low. (2) The slope of the diagonal patterns of the N 0 time-height field is low (because Atmos. Meas. Tech., 9,[3145][3146][3147][3148][3149][3150][3151][3152][3153][3154][3155][3156][3157][3158][3159][3160][3161][3162][3163]2016 www.atmos-meas-tech.net/9/3145/2016/ N 0 is driven by small drops, with low terminal velocities).
(3) Patches of raindrops are formed by upward winds (see Fig. 9), and these patches disappear as soon as the wind in the lower layers turns downward. For instance, around 19:55, there is an upward wind between 100 and 400 m, which results in a patch of drops between 300 and 500 m. We can also note that even if the absolute values of N 0 are dubious, these patches are probably real, since they are driven by the wind. Because the DSD for this event is globally exponential (see Fig. 8), it is reasonable to estimate a DSD slope (Fig. 9, bottom). In addition, because this slope is calculated over the entire range of diameters, the retrieved slope field is vertically consistent (less affected by the wind). We also note that the retrieved slopes are consistent with the results of Giangrande et al. (2010) (Fig. 8 in this paper) in a similar case study.

Conclusions
We have built a 4-D-VAR data assimilation algorithm for retrieving vertical DSD profiles and vertical wind under the bright band from observations coming from a micro-rain radar (MRR) and from a co-located disdrometer, associated with a vertical propagation model. In this paper, we focused on the data assimilation algorithm. The algorithm finely handles measurements of various natures collected at different resolutions. We consequently chose to use a simple propagation model, only taking gravity and vertical wind into account. Because of the limitations of the model, the retrieval algorithm is currently only suitable to study stratiform, light rain events.
The algorithm was firstly applied to simulated data in a twin experiment context. In simulated observations, despite the addition of instrumental noise and turbulence effects, we showed that the proposed technique was able to retrieve DSD fields along with vertical wind patterns and intensities thanks to radar Doppler spectra and disdrometer drop fluxes. In particular, the algorithm appears able to handle the turbulence impact on Doppler spectra even if the turbulence is not explicitly implemented in the propagation model.
The behavior of the algorithm is assessed on a light rain event during which co-localized MRR and optical disdrometer were used in the south of France during the HyMeX campaign. According to this work, the algorithm appears able to make the best use of the data from the two instruments. In particular, we noticed (suspicious) vertical trends on the MRR products which cannot be reproduced by the propagation model due to its conservative nature. However, despite this trend, the assimilation algorithm is able to produce a good compromise between the observations at different heights. The vertical retrieved wind field, although not independently validated, is spatially and temporally coherent with what can be expected for a rain event such as the one studied here. The DSD fields also appear spatially coherent and self-consistent. By combining ground and radar data, the presented algorithm is therefore able to retrieve the vertical wind field and to improve the DSD profile retrieval.
Independent information on the vertical wind and/or on the DSD are nevertheless necessary to estimate the performances of the algorithm. For this purpose, we plan to use dual-frequency wind profilers (see for instance Williams, 2012, which retrieves vertical wind in rain from two profilers). Another possibility would be to use the methods developed by Giangrande et al. (2010) or Tridon and Battaglia (2015), which retrieve vertical winds and DSD from K band and/or W band Doppler spectra. We could also study the impact on the quality of the retrievals of the addition of other frequencies directly in the assimilation process.
Some improvements are needed to provide an algorithm suitable for various weather situations (tropical rain and convective events, for instance). In our model, we discarded horizontal motions. This hypothesis is valid only if there is no significant horizontal wind or at least if the DSD is sufficiently homogeneous in the horizontal plane. It will be necessary to evaluate the impact of such an hypothesis more precisely. Another possibility would consist in taking into account the horizontal wind by extending the model with an additional spatial dimension. In our case, this can be done with the use of the YAO assimilation tool that provides a simple way to do it.
Other processes such as evaporation and coalescence/breakup have to be implemented in the propagation model. This enhanced model would allow the relative importance of the different physical phenomena (wind, evaporation, collisions, . . . ) to be investigated, and therefore which are really essential for the assimilation algorithm to be determined. It would also help to better parameterize the algorithm through its ability to explain the observations. Finally, the method is well suited for merging observations of different types each with its own temporal and spatial resolution. Indeed, when merging new observations, only the observation operator has to be adapted. For instance we could use more disdrometers located in the vicinity or the reflectivity at another frequency from another radar. Any kind of sensor which observes geophysical parameters playing a role during the fall of raindrops could also be easily introduced in the assimilation algorithm (ground wind, pressure, humidity).

Data availability
The data sets used in this paper (for instance micro-rain radar measurements) are available on the HyMeX website: http://mistrals.sedoo.fr/HyMeX/.