3D trajectories and velocities of rainfall drops in a multifractal turbulent wind field

Weather radars measure rainfall in altitude, whereas hydro-meteorologists are mainly interested in rainfall at ground level. During their fall, drops are advected by the wind, which affects the location of the measured field.

The governing equation of a rain drop's motion relates the acceleration to the forces of gravity and buoyancy along with the drag force. It depends non-linearly on the instantaneous relative velocity between the drop and the local wind, which yields complex behaviour. Here, the drag force is expressed in a standard way with the help of a drag coefficient expressed as a function of the Reynolds number. Corrections accounting for the oblateness of drops greater than 1–2 mm are suggested and validated through a comparison of the retrieved “terminal fall velocity” (i.e. without wind) with commonly used relationships in the literature.

An explicit numerical scheme is then implemented to solve this equation for a 3+1D turbulent wind field, and hence analyse the temporal evolution of the velocities and trajectories of rain drops during their fall. It appears that multifractal features of the input wind are simply transferred to the drop velocity with an additional fractional integration whose level depends on the drop size, and a slight time shift. Using an actual high-resolution 3D sonic anemometer and a scale invariant approach to simulate realistic fluctuations of wind in space, trajectories of drops of various sizes falling form 1500 m are studied. For a strong wind event, drops located within a radar gate in altitude during 5 min are spread on the ground over an area of the size of a few kilometres. The spread for drops of a given diameter is found to cover a few radar pixels. Consequences on measurements of hydro-meteorological extremes that are needed to improve the resilience of urban areas are discussed.

1 Introduction

During their fall, drops are advected by wind. Quantitative rainfall estimation with the help of weather radars is affected by this issue since drops can be displaced horizontally between their measurement location in altitude and their ground impact location, which is of interest for hydro-meteorologists. This effect is usually called wind drift in the literature and sometimes wind advection. The potential bias and uncertainty introduced in radar measurements is stronger at higher resolution, i.e. typically with pixel sizes smaller than 1–2 km², which are needed for urban applications, for example. Collier (1999) suggests that correction schemes should be implemented for this kind or higher radar resolution. Lauri et al. (2012) reported that far from the radar (i.e. typically more than 150 km), even with low elevation (0.3^∘), displacements of few tens of kilometres are found, which actually distort the measured area.

Most correction schemes rely on the use of 4D wind profiles derived from numerical prediction models (Mittermaier et al., 2004; Lack and Fox, 2007; Lauri et al., 2012; Sandford, 2015) or a combination of them with reanalysis (Dai et al., 2013, 2019; Yang et al., 2020). The latter also accounts for drop size distribution (DSD). They report an improvement by ≈3 % of the correlation between radar and rain gauge measurement and a reduction of discrepancy of ≈18 % over eight selected events. Lack and Fox (2007) directly used Doppler radar wind measurement at 2.5 km scale to adjust for the wind drift effect. In general, correction schemes use wind data of rather coarse resolution (typically km(s)) and assume a constant wind shear. Nevertheless, some variability at smaller space-time scales is usually acknowledged, especially during convective events, i.e. those for which wind drift causes the greatest uncertainty (Lack and Fox, 2007).

Wind effects on rainfall drops is also reported to generate discrepancies between the vertical velocities measured and expected terminal fall velocities. For example, Montero-Martínez and García-García (2016) studied events with calm, light and moderate wind with various rainfall levels and found a widening of the fall velocity distribution under windy conditions. They found super-terminal drops only for diameters <0.7 mm and more often under wind conditions. Sub-terminal fall velocities for drops of sizes up to 2 mm are reported. Bringi et al. (2018) found that under low wind speed and turbulence, no discrepancies with expectations are found, while under high wind speed and turbulence, there is a clear widening of the distribution. A linear decrease of the mean fall velocity with increasing turbulence intensity is reported. Maximum decreases of 25 %–30 % are observed. Thurai et al. (2019) also found such decreases for drops greater than 2 mm in high turbulence intensity conditions. It is associated to an asymmetry that also appears in the drop shape. They also found that horizontal drop velocities in both direction and magnitude show a “remarkable agreement” with the wind sensor at 10 m. Stout et al. (1995) explored the effect of the nonlinear drag coefficient on the fall velocity through numerical simulations. They showed that even heavy drops exhibited a reduced settling velocity in isotropic turbulence.

Turbulence is found to have contradictory effects on the distribution of the fall velocity. Indeed, increasing the turbulence level in windy and rainfall conditions will yield more collision and breakup, resulting in smaller drops inheriting the speed of larger parent drops, and hence observations of super-terminal velocities. On the other hand, turbulence is said to yield a decrease in fall velocities because drops (especially ones of <1 mm) are more affected by eddies.

Such findings on the discrepancies between observed and expected fall velocities have effects on the relation between rainfall and kinetic energy, i.e. the erosivity “power” of rainfall (Pedersen and Hasholt, 1995) and also building performance to outdoor conditions (Tian et al., 2018; Blocken et al., 2011).

The studies previously mentioned basically do not account for small-scale wind fluctuations in both space and time. In this paper, we suggest studying the behaviour of individual rainfall drops of various sizes in a high-resolution turbulent wind field. The variability of the wind is accounted for through the framework of universal multifractals (UM) (see Schertzer and Tchiguirinskaia, 2020, for a recent review). Such a physically based framework is designed to analyse and simulate geophysical fields exhibiting extreme variability over a wide range of space-time scales as wind. Drop oblateness is also accounted for.

The paper is organised as follows. In Sect. 2, a deterministic equation for the fall of individual oblate drops in a 3D field is derived and validated through the comparison of the terminal fall velocity obtained with commonly used formulas. In Sect. 3, the framework of universal multifractals is described briefly. Then, the drops are subjected to simulated multifractal fields as the wind input and multifractal behaviour of the horizontal drop velocity is assessed. Finally, in Sect. 4, 3D wind is reconstructed from high-resolution 3D sonic anemometer data and strong scaling assumptions. This field is used to study the trajectories of drops between falling from 1500 m to the ground.

2 A deterministic equation for oblate drops in a wind field

2.1 Formulation of the equation

Let us denote $(x,y,z)$ the horizontal, lateral and vertical coordinates in a standard Cartesian framework with unit vectors (${\underset{\mathrm{¯}}{e}}_x$, ${\underset{\mathrm{¯}}{e}}_y$, ${\underset{\mathrm{¯}}{e}}_z$). We aim at writing the motion equation of a particle of water (a drop) of velocity ${\underset{\mathrm{¯}}{v}}_p$, density ρ_p and falling in the atmosphere under the influence of the gravity $\underset{\mathrm{¯}}{g}=-g{\underset{\mathrm{¯}}{e}}_z$ (where g=9.81 m s⁻²) and the wind ${\underset{\mathrm{¯}}{v}}_{\mathrm{wind}}$ (with three components). The density of the atmosphere is denoted ρ_air. The water particle is characterised by its equivolumic diameter D_eq, which corresponds to the diameter of the sphere having the same total volume. Hence we have Vol $=\frac{\mathit{\pi }}{\mathrm{6}}{D}_{\mathrm{eq}}^{\mathrm{3}}$. Finally, the relative velocity between the wind and the falling particle is ${\underset{\mathrm{¯}}{v}}_{\mathrm{rel}}={\underset{\mathrm{¯}}{v}}_{\mathrm{wind}}-{\underset{\mathrm{¯}}{v}}_p$ (${\underset{\mathrm{¯}}{v}}_{\mathrm{rel}}$ is the vector with three components, while v_rel is its Euclidean norm).

The drop is subjected to three forces:

• The gravity equal to ${\mathit{\rho }}_p\text{Vol}\underset{\mathrm{¯}}{g}$.

• The buoyancy equal to $-{\mathit{\rho }}_{\mathrm{air}}\text{Vol}\underset{\mathrm{¯}}{g}$.

• The drag, which is commonly written as $\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathit{\pi }{D}^{\mathrm{2}}}{\mathrm{4}}{c}_{\mathrm{D}}{\mathit{\rho }}_{\mathrm{air}}{v}_{\mathrm{rel}}{\underset{\mathrm{¯}}{v}}_{\mathrm{rel}}$. Re is the common Reynolds number $\mathit{Re}=\frac{{\mathit{\rho }}_{\mathrm{air}}{v}_{\mathrm{rel}}D}{{\mathit{\mu }}_{\mathrm{air}}}$, where μ_air is the absolute viscosity of air; c_D is the drag coefficient and depends in general on Re and D_eq. The next section is devoted to its determination.

As a consequence, the equation of motion of the falling particle is given by Newton's second law, which equals the mass times the acceleration to the net force (here, it was divided by the mass):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}{\underset{\mathrm{¯}}{v}}_p}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{3}}{\mathrm{4}\text{D}}}{c}_{\mathrm{D}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{air}}}{{\mathit{\rho }}_p}}{v}_{\mathrm{rel}}{\underset{\mathrm{¯}}{v}}_{\mathrm{rel}}+\underset{\mathrm{¯}}{g}{\displaystyle \frac{{\mathit{\rho }}_p-{\mathit{\rho }}_{\mathrm{air}}}{{\mathit{\rho }}_p}}.\end{array}$$

2.2 Determination of the drag coefficient

Before discussing how the drag coefficient is determined, it should be recalled that the rainfall drops considered in this paper are not spherical. Indeed, drops greater than typically 1.5 mm become oblate in their fall. This oblateness increases with size. A very commonly used model consists in an ellipsoid with an axis ratio varying depending on the size. Thurai et al. (2007) showed that such model is too simplistic since drops are not symmetric in the direction perpendicular to their fall. Following an in-depth analysis of the drop shape assessed with the help of a 2D-video disdrometer (Kruger and Krajewski, 2002) in the measurement campaign of an artificial rainfall experiment, they suggested the following formula for the shape:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}x& =a_{\mathrm{1}}\sqrt{\mathrm{1}-{\left({\displaystyle \frac{z}{a_{\mathrm{2}}}}\right)}^{\mathrm{2}}}\left[{\cos }^{-\mathrm{1}}\left({\displaystyle \frac{z}{a_{\mathrm{2}}{a}_{\mathrm{3}}}}\right)\right]\left[a_{\mathrm{4}}{\left({\displaystyle \frac{z}{a_{\mathrm{2}}}}\right)}^{\mathrm{2}}+\mathrm{1}\right]\\ & =f\left(z\right)\end{array}\end{array}$$

with:

$$\begin{array}{}\text{(3)}& \begin{array}{lc}{a}_{\mathrm{1}}=& \frac{\mathrm{1}}{\mathit{\pi }}\left(\mathrm{0.02914}{D}_{\mathrm{eq}}^{\mathrm{2}}+\mathrm{0.9263}{D}_{\mathrm{eq}}+\mathrm{0.07791}\right)\\ a_{\mathrm{2}}=& -\mathrm{0.01938}{D}_{\mathrm{eq}}^{\mathrm{2}}+\mathrm{0.4698}{D}_{\mathrm{eq}}+\mathrm{0.09538}\\ a_{\mathrm{3}}=& -\mathrm{0.06123}{D}_{\mathrm{eq}}^{\mathrm{3}}+\mathrm{1.3880}{D}_{\mathrm{eq}}^{\mathrm{2}}-\mathrm{10.41}{D}_{\mathrm{eq}}^{\mathrm{2}}+\mathrm{28.34}\\ a_{\mathrm{4}}=& -\mathrm{0.01352}{D}_{\mathrm{eq}}^{\mathrm{3}}+\mathrm{0.2014}{D}_{\mathrm{eq}}^{\mathrm{2}}-\mathrm{0.8964}{D}_{\mathrm{eq}}^{\mathrm{2}}\\ & +\mathrm{1.226}\text{if}{D}_{\mathrm{eq}}>\mathrm{4}\mathrm{mm}\\ a_{\mathrm{4}}=& \mathrm{0}\text{if}\mathrm{1.5}\mathrm{mm}\le D_{\mathrm{eq}}\le \mathrm{4}\mathrm{mm}.\end{array}\end{array}$$

This shape corresponding to a solid of revolution around the z axis is used in this paper. It is displayed in Fig. 2a for drops with equivolumic diameter ranging from 1.5 to 5.5 mm. It should be mentioned that computing the volume as an integral of the shape (Vol $={\int }_{z_{\min }}^{z_{\max }}\mathit{\pi }f(z{)}^{\mathrm{2}}\mathrm{d}z$, see Fig. 1) yields minor differences with the expected volume of $\frac{\mathit{\pi }{D}_{\mathrm{eq}}^{\mathrm{3}}}{\mathrm{6}}$. They are highlighted in Fig. 2b. In order to account for this small difference, once an equivolumic diameter is set, the corresponding one that would lead to the expected volume from Eqs. (2) and (3) is computed from a correspondence table. The relationship, which is obviously close to the bisector, is displayed in Fig. 2c with the horizontal axis corresponding to the real D_eq of the drop and the vertical axis to the D_eq to be input in Eqs. (2) and (3) to retrieve the correct expected volume. A consequence is that the oblateness of drops will be considered only from an equivolumic diameter greater than 1.527 mm.

Figure 1 Illustration of how the volume (Vol), the surface area (SA) and the mean projected longitudinal cross-sectional area (MPA_∥) of the particle used in this paper can be computed via integral calculus. The drop is actually a solid of revolution around the “z” axis.

Figure 2 (a) Drop shape model used in this paper; (b) drop volume vs. equivolumic diameter; (c) diameter relation to retrieved wanted volume; (d) parameters characterising non-spherical shape of drops vs. equivolumic diameter; (e) drag coefficient c_D vs. Re number and (f) terminal fall velocity vs. equivolumic diameter.

For non-spherical shapes, it is quite tricky to compute the corresponding drag coefficient as a function of the Reynolds number. The literature about this issue is quite abundant, and the interested reader is referred to Chap. 4 of the PhD dissertation of Bagheri (2015) or Hölzer and Sommerfeld (2008) for details. In the approach that they implemented, three parameters are used to characterize the non-spherical shapes of the falling particle with the help of three dimensionless parameters: the sphericity, the crosswise sphericity and the lengthwise sphericity. The latter two depend on the orientation of the particle with regards to the flow. Here, it is assumed that drops are oriented perpendicularly to the flow, i.e. the “z” axis of Eq. (2) is parallel to ${\underset{\mathrm{¯}}{v}}_{\mathrm{rel}}$. In the general case, these parameters may be complex to assess but with the shape derived from Eq. (2) (Thurai et al., 2007), theoretical formula can be obtained. See Fig. 1 for an illustration of the computations via integral calculus. The three parameters are:

• The sphericity ψ, which is equal to ratio between the surface area (SA) of the equivolumic sphere to the actual surface area of the particle. It is equal to 1 for sphere and decreases for increasingly and fewer spherical particles; $\mathit{\psi }=\frac{pi{D}_{\mathrm{eq}}^{\mathrm{2}}}{\text{Surface Area}}$. That is, within the framework of this paper, we have SA $={\int }_{z_{\min }}^{z_{\max }}\mathrm{2}\mathit{\pi }f\left(z\right)\sqrt{\mathrm{1}+f^{\prime }(z{)}^{\mathrm{2}}}\mathrm{d}z$.

• The crosswise sphericity ψ_⟂, which is equal to the ratio between the projected area of the volume equivalent sphere and the projected area of the particle normal to the falling direction (here, ${\underset{\mathrm{¯}}{e}}_z$); ${\mathit{\psi }}_{⟂}=\frac{D_{\mathrm{eq}}^{\mathrm{2}}}{D_{\max }^{\mathrm{2}}}$. It is equal to 1 for sphere and decreases for larger drops since they become oblate.

• The lengthwise sphericity ψ_∥, which is defined as the cross-sectional area of the volume equivalent sphere divided by the difference between half the surface area and the mean projected longitudinal cross-sectional area of particle (MPA _∥); ${\mathit{\psi }}_{\parallel }=\frac{\mathit{\pi }{\left(\frac{D_{\mathrm{eq}}}{\mathrm{2}}\right)}^{\mathrm{2}}}{\frac{\text{SA}}{\mathrm{2}}-{\text{MPA}}_{\parallel }}$. In the specific drop model of this paper, we have MPA ${}_{\parallel }={\int }_{z_{\min }}^{z_{\max }}\mathrm{2}f\left(z\right)\mathrm{d}z$ (it basically corresponds to the 2D area of the drop plotted in Fig. 1).

The evolution of these parameters as a function of D_eq for the drops considered is shown in Fig. 2d. The increasing oblateness of drops with increasing size is translated through the fact that the parameters are getting further away from 1. In order to define the drag coefficient, the corrections suggested by Hölzer and Sommerfeld (2008) to account for non-sphericity of particles are then implemented in the formula of White (1974), previously used by Stout et al. (1995) who worked only on spherical drops. This yields

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{c}_{\mathrm{D}}=& {\displaystyle \frac{\mathrm{8}}{\mathit{Re}\sqrt{{\mathit{\psi }}_{\parallel }}}}+{\displaystyle \frac{\mathrm{16}}{\mathit{Re}\sqrt{{\mathit{\psi }}_{⟂}}}}+{\displaystyle \frac{\mathrm{6}}{\left(\mathrm{1}+\sqrt{\mathit{Re}}\right){\mathit{\psi }}^{\mathrm{3}/\mathrm{4}}}}\\ & +{\displaystyle \frac{\mathrm{0.25}\times {\mathrm{10}}^{\mathrm{0.4}(-\log \mathit{\psi }{)}^{\mathrm{0.2}}}}{{\mathit{\psi }}_{⟂}}}.\end{array}\end{array}$$

The evolution of c_D as a function of Re for various drop parameters is displayed in Fig. 2e and follows standard patterns.

2.3 Validation of the formula

In order to validate the developed equation, the retrieved terminal fall velocity is assessed for each equivolumic diameter. It corresponds to the velocity of the permanent regime with no wind, i.e. the drag plus the buoyancy exactly compensate the gravity. Computations are carried out with ρ_air=1.205 kg m⁻³, ${\mathit{\mu }}_{\mathrm{air}}=\mathrm{1.81}\times {\mathrm{10}}^{-\mathrm{5}}$ kg m⁻¹ m⁻², ρ_water=998.2 kg m⁻³ g=9.81 m s⁻², as in Stout et al. (1995).

The relation between the terminal fall velocity obtained vs. the equivolumic diameter is displayed in red in Fig. 2f. The developed equations enable us to retrieve the commonly used relation (Beard, 1977; Lhermitte, 1988; Best, 1950; Atlas et al., 1973) for drops of diameters of up to 4 mm. The deviations found when considering spherical drops (in green) are visible for diameters greater than 2 mm, which highlights the need to account for drop oblateness.

2.4 Numerical scheme for solving the equation

Equation (1) is solved numerically through the implementation of a simple Eulerian numerical scheme. Within such a framework: (i) a discretization of time with time step Δt is introduced yielding discrete time steps $t_n=n\times \mathrm{\Delta }t$, where n is an integer; (ii) we aim at finding an approximation of ${\underset{\mathrm{¯}}{v}}_p$ at time step n denoted ${\underset{\mathrm{¯}}{v}}_{p,n}$; (iii) the first derivative in Eq. (1) is approximated as $\frac{\mathrm{d}{\underset{\mathrm{¯}}{v}}_p}{\mathrm{d}t}\left(t_n\right)\approx \frac{{\underset{\mathrm{¯}}{v}}_{p,n+\mathrm{1}}-{\underset{\mathrm{¯}}{v}}_{p,n}}{\mathrm{\Delta }t}$. This leads to the following equation for the numerical scheme:

$$\begin{array}{}\text{(5)}& {\underset{\mathrm{¯}}{v}}_{p,n+\mathrm{1}}={\underset{\mathrm{¯}}{v}}_{p,n}+\mathrm{\Delta }t\left[{\displaystyle \frac{\mathrm{3}}{\mathrm{4}D}}{c}_{\mathrm{D},n}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{air}}}{{\mathit{\rho }}_p}}{v}_{\mathrm{rel},n}{\underset{\mathrm{¯}}{v}}_{\mathrm{rel},n}+\underset{\mathrm{¯}}{g}{\displaystyle \frac{{\mathit{\rho }}_p-{\mathit{\rho }}_{\mathrm{air}}}{{\mathit{\rho }}_p}}\right],\end{array}$$

where ${\underset{\mathrm{¯}}{v}}_{\mathrm{rel},n}$ and c_D,n are computed at time step t_n using the formulas discussed in Sects. 2.1 and 2.2. Assuming some initial conditions (always no horizontal velocity and a vertical one equal to the terminal fall for the corresponding diameter), it is then possible to reconstruct the time series of velocity for the drops. From this, the temporal evolution of the position (i.e. the trajectory) is derived. It is necessary to properly assess the wind accounting for the current position of the drop. A time step of Δt=0.01 s is used in this paper, and it was checked that it ensured a stability of the numerical scheme.

3 Behaviour of horizontal drop velocity with multifractal input

3.1 Brief recollection of the universal multifractal framework

It is outside the scope of the paper to introduce the framework of universal multifractals (UM) in detail. Hence, only the most important elements are recalled here, and interested readers are referred to the references mentioned or to a recent review by Schertzer and Tchiguirinskaia Schertzer and Tchiguirinskaia (2020) for more details.

Let us consider a field ϵ_λ at a resolution λ defined as the ratio between the outer scale (L) and the observation scale (l), $\mathit{\lambda }=L/l$. For multifractal fields, the moment of order q of the field is a power law related to the resolution:

$$\begin{array}{}\text{(6)}& \langle {\mathit{ϵ}}_{\mathit{\lambda }}^q\rangle \approx {\mathit{\lambda }}^{K\left(q\right)},\end{array}$$

where K(q) is the scaling moment function. It fully characterizes the variability across scales of the field. Within the specific framework of UM (Schertzer and Lovejoy, 1987, 1997), towards which multiplicative cascades processes converge, only two parameters with physical interpretation are needed to characterize K(q) for conservative fields:

• C₁, the mean intermittency co-dimension, which measures the clustering of the (average) intensity at smaller and smaller scales; C₁=0 for an homogeneous field;

• α, the multifractality index ($\mathrm{0}\le \mathit{\alpha }\le \mathrm{2}$), which measures the clustering variability with regards to the intensity level.

For UM, we have:

$$\begin{array}{}\text{(7)}& K\left(q\right)={\displaystyle \frac{C_{\mathrm{1}}}{\mathit{\alpha }-\mathrm{1}}}(q^{\mathit{\alpha }}-q).\end{array}$$

A non-conservative field (ψ_λ), i.e. whose mean is not preserved across the scale can be written as ${\mathit{\psi }}_{\mathit{\lambda }}\approx {\mathit{ϵ}}_{\mathit{\lambda }}{\mathit{\lambda }}^{-H}$, where H is the non-conservativeness parameter; H=0 for conservative fields. Positive values correspond to a fractional integration to go from ϵ_λ to ψ_λ and to stronger correlations within the field ψ_λ. Negative values correspond to a fractional differentiation; H is typically between 0 and 1 for geophysical fields.

The first step of a multifractal analysis usually consists in a spectral analysis. For multifractal fields, the power spectra (E) should scale with wave number k:

$$\begin{array}{}\text{(8)}& E\left(k\right)=k^{-\mathit{\beta }},\end{array}$$

with the spectral slope β

$$\begin{array}{}\text{(9)}& \mathit{\beta }=\mathrm{1}+\mathrm{2}H-K_c\left(\mathrm{2}\right),\end{array}$$

where K_c is the scaling moment function (Eq. 7) of the conservative part of the field. To analyse the latter, a trace moment (TM) is implemented. It notably enables us to assess the quality of the scaling behaviour. It basically consists in plotting Eq. (6) in log–log. Straight lines should be retrieved, and the slope gives K(q). Finally, UM parameters are estimated with the help of the double trace moment (DTM) technique, which is tailored for UM fields and enables robust estimation of UM parameters (Lavallée et al., 1993).

3.2 Methodology

In this section, the scaling behaviour of the horizontal drop velocity is assessed using numerical simulations. Working with such input whose features are fully known is helpful to understand how drops react to wind.

More precisely, a horizontal input v_x,wind for Eq. (1) is simulated with the help of blunt multifractal discrete cascades (Gires et al., 2020). Such a process yields only positive values, which is not realistic for wind. Hence, a standard “complex trick” was used to generate a field with both positive and negative values (Schertzer and Lovejoy, 1995). To implement it, two fields X₁ and X₂ are generated with the wanted features, and a third one is obtained with the help of the following equation (Real is the real part):

$$\begin{array}{}\text{(10)}& X=\text{Real}\left[\exp \left(\log X_{\mathrm{1}}+i\log X_{\mathrm{2}}\right)\right].\end{array}$$

Such a field divided by 2 was used as input. With UM parameter α=1.7 and C₁=0.2 1024 long time step series are generated, which corresponds to typical values for turbulent wind fields (Fitton et al., 2011). The time step is assumed to be 0.01 s, which means that drops are basically studied over 10 s. For the initial conditions, drops are assumed to have no horizontal velocity and a vertical component equal to its corresponding terminal fall velocity. Since scaling is a statistical behaviour, an ensemble of 100 independent samples was generated, and the corresponding ensemble of horizontal drop velocity was simulated using Eq. (5) for drops of various sizes ($D_{\mathrm{eq}}\subset [\mathrm{0.1},\mathrm{0.2},\mathrm{0.4},\mathrm{0.6},\mathrm{0.8},\mathrm{1},\mathrm{1.5},\mathrm{2},\mathrm{3},\mathrm{4}]$).

3.3 Results and discussion

Figure 3 displays the temporal evolution of drops' horizontal velocity over 10 s for a sample of wind input (shown in black). Three drop diameters are displayed (0.1, 0.6 and 2 mm). It can be seen, notably on the zoomed in part of the figure (lower panel) that the smaller drop (D_eq=0.1 mm, shown in blue) follows wind fluctuations well, with only a limited dampening of the fluctuations. A small delay (≈0.01 s) corresponding to a reaction time is noted. As can be expected, larger drops (D_eq=0.6 mm, shown in green; and D_eq=2 mm, shown in red) tend to dampen wind fluctuations even more.

Figure 3 (a) Temporal evolution of the drop horizontal velocity during 10 s for various drop diameters (the equivolumic diameter is indicated) with the same multifractal wind input (black). (b) A zoom of the above curve on a shorter period.

In order to quantify this qualitative behaviour more precisely, a multifractal analysis on the retrieved ensembles was performed. Figure 4 displays the outcome of spectral and TM analyses for drops of equivolumic diameter ranging from 0.1 to 2 mm. The spectral analyses reflect a good scaling behaviour over the whole range of scales. Spectral slopes (β in Eq. 8) of 0.86 and 2.25, respectively, are retrieved. For the 2 mm drop, the value corresponds to non-conservative fields. In order to ensure that a conservative field is studied in TM analysis, which is necessary (Lavallée et al., 1993), a fractional differentiation with an exponent $(\mathit{\beta }-\mathrm{1})/\mathrm{2}$ is implemented on the field before implementing this TM analysis. TM analysis is displayed in the right-hand column of Fig. 4. For the 0.1 mm drop, an excellent scaling behaviour is retrieved with the coefficient of determination r² for q=1.5 greater than 0.99. DTM analysis yields estimates of UM parameters α, C₁ and H equal to 1.68, 0.21 and 0.12, respectively, which is close to the features of the input series. For the 2 mm drop, the scaling is slightly degraded but remains good (r²=0.95 for q=1.5), and we find α=1.69, C₁=0.14 and H=0.79.

Figure 4 Scaling behaviour of the simulated drop velocity for D_eq equal to 0.1 mm (top row) and 2 mm (bottom row). Spectral analysis (i.e. Eq. 8 in log–log) on the direct simulations is displayed in the left-hand column. TM analysis (i.e. Eq. 6 in log–log) on the velocity after implementing a fractional integration is displayed in the right-hand column.

Figure 5 displays a summary of the UM analysis carried out on the generated series for the various drops. The scaling behaviour is excellent for small drops and remains good for all drop sizes with r² for q=1.5 always greater than 0.95 (Fig. 5e). The need for a fractional differentiation before implementing TM analysis is visible with the very poor scaling found when analysing the field directly. The non-conservativeness parameter rapidly increases from 0.1 to 0.8 with the drop size increasing from 0.1 to ≈ 1–1.5 mm. For larger drops, it remains rather stable. This increase of H is basically a quantification of the increased dampening of wind fluctuations observed for larger drops, shown in Fig. 3. With respect to the UM parameters α and C₁, the former remains stable and close to the input value of 1.7 for all drop sizes. The latter exhibits a small decrease with larger drops. It should be recalled that this approach is somehow artificial since all drops perceive the same wind, which would not be the case in reality because they do not fall at the same vertical speed. In summary, this investigation shows that horizontal drop velocity basically reproduces the multifractal properties of the wind input with an increased level of non-conservativeness H; H increases strongly for drops smaller than 1 mm and then stabilizes.

Figure 5 Summary of the multifractal analysis performed on the ensembles of simulated horizontal drop velocities using a wind input with α=1.7, C₁=0.2. The various multifractal parameters are displayed vs. D (the equivolumic drop diameter).

4 Ground impact location of drops falling in a turbulent wind field

4.1 Methodology

The purpose of this section is to investigate where drops falling from a height of 1500 m reach the ground. Given the time step of 0.01 s used in the equation and the fact that drops move in space during their fall, this means that it is necessary to have high-resolution space-time 3D wind data over an area of the typical size of a few kilometres to fully address the issue. Such data is unfortunately not available. Hence, we suggest here to reconstruct a somehow realistic wind from a punctual measurement relying on previous findings on turbulence.

4.1.1 Wind data

More precisely, we use 100 Hz 3D sonic anemometer data collected at by a device installed at 78 m on a meteorological mast located on the Pays d'Othe wind farm within the framework of the ANR RW-Turb project. The wind farm is roughly 120 km south-east of Paris on a slightly sloppy area. More details can be found in the data paper under discussion at ESSD (Gires et al., 2022) and in the data set (Gires et al., 2021). Two wind series with very different average horizontal wind speeds (i.e. 1.8 vs. 11.8 m s⁻¹) were extracted for this study. They were later denoted low wind event and strong wind event. The corresponding time series of the three components of the wind for both events are displayed in Fig. 6 over approximately 900 s. The low wind event was collected on 20 January 2021, while the strong one occurred on 6 January 2021.

Figure 6 Temporal evolution of the 100 Hz wind data from a 3D sonic anemometer for the low (top) and strong (bottom) wind events used in this paper.

4.1.2 Generation of an anisotropic 3D turbulent field

In this section, we discuss how to stochastically generate a turbulent field reproducing the physics of such a flow constraint as well as possible to have the empirical velocity values $v(x,y,z,t)$ of the 3D sonic anemometer located at $(x,y,z,t)$.

Scaling and anisotropy

It is quite obvious that gravity has such a strong impact on drop trajectories and dynamics of drops that classical scaling approaches just fail because they presuppose isotropy. On the contrary, the anisotropy between the vertical and the horizontal induced by gravity is so ubiquitous in geophysics that it has led to the general concept and framework of “generalized scale invariance” for the analysis and simulation of anisotropic fields (Schertzer and Lovejoy, 1985, 1988, 1989; Lazarev et al., 1994; Schertzer et al., 2012), and new extensions have been developed for vector fields (Schertzer and Tchiguirinskaia, 2015, 2020). While classical approaches consider scaling only after assuming isotropy, general scale invariance first posit scaling and then study the remaining non-trivial symmetries. Because this paper deals with scalar discrete cascades, and, therefore, their limitations, we consider an oversimplification of the generalised scale invariance. Recall that for the simplest case of generalised scale invariance, a v horizontal component of a statistically translation invariant velocity field satisfies the following self-affine scale invariance:

$$\begin{array}{}\text{(11)}& \mathrm{\Delta }v\left(T_{\mathit{\lambda }}\mathrm{\Delta }\underset{\mathrm{¯}}{x}\right)\stackrel{d}{=}{\mathit{\lambda }}^{H_{\mathrm{h}}}\mathrm{\Delta }v\left(\mathrm{\Delta }\underset{\mathrm{¯}}{x}\right),\end{array}$$

where $\stackrel{d}{=}$ stands for equality in distribution, $\mathrm{\Delta }\underset{\mathrm{¯}}{x}$ for a (vector) pair separation, $\mathrm{\Delta }v=\left|v(\underset{\mathrm{¯}}{x}+\mathrm{\Delta }\underset{\mathrm{¯}}{x})-v\left(\underset{\mathrm{¯}}{x}\right)\right|$ for the induced velocity component shear, H_h for the horizontal scaling exponent ($H_{\mathrm{h}}=\mathrm{1}/\mathrm{3}$ for Kolmogorov's scaling) and T_λ for a self-affine change of scale:

$$\begin{array}{}\text{(12)}& T_{\mathit{\lambda }}={\mathit{\lambda }}^G=\exp \left(\log \right(\mathit{\lambda }G);G=\left(\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& H_z\end{array}\right),\end{array}$$

where G is the generator of T_λ, $H_z=H_{\mathrm{h}}/H_{\mathrm{v}}$ is the scaling anisotropy exponent and H_v the scaling exponent along the vertical; $H_{\mathrm{v}}=\mathrm{3}/\mathrm{5}$ for a Bolgiano–Obukhov scaling along the vertical, and, therefore, $H_z=\mathrm{5}/\mathrm{9}$, while for isotropic 3D turbulence $H_z=\mathrm{1},H_{\mathrm{h}}=H_{\mathrm{v}}$. The main property of T_λ is that it broadly generalises in a straightforward manner scales $|.|$ from classical isotropic frameworks (H_z=1) to generalised scales $\Vert .\Vert$ of anisotropic frameworks (e.g. H_z≠1), so that any isotropic cascade can be transformed in this manner into an anisotropic cascade without anything else. Just recall the canonical example of generalised scales for a diagonal generator (p≥1):

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{G}_{i,j}=& g_i{\mathit{\delta }}_{i,j}:\Vert \underset{\mathrm{¯}}{x}{\Vert }_p=\left(\sum _i\right|x_i{]}^{p/g_i}{)}^{\mathrm{1}/p}\\ & \Rightarrow \Vert T_{\mathit{\lambda }}\underset{\mathrm{¯}}{x}{\Vert }_p=\mathit{\lambda }\Vert \underset{\mathrm{¯}}{x}{\Vert }_p,\end{array}\end{array}$$

which displays the (generalised) scaling property of the generalised scales.

Equation (11) is valid for any direction of $\mathrm{\Delta }\underset{\mathrm{¯}}{x}$ ($\underset{\mathrm{¯}}{x}=(x,y,z)$), in particular along the horizontal Δx and the vertical Δz:

$$\begin{array}{}\text{(14)}& \mathrm{\Delta }v\left(\mathit{\lambda }\mathrm{\Delta }x\right)\stackrel{d}{=}{\mathit{\lambda }}^{H_{\mathrm{h}}}\mathrm{\Delta }v\left(\mathrm{\Delta }x\right);\mathrm{\Delta }v\left(\mathit{\lambda }\mathrm{\Delta }z\right)\stackrel{d}{=}{\mathit{\lambda }}^{H_{\mathrm{v}}}\mathrm{\Delta }v\left(\mathrm{\Delta }v\right),\end{array}$$

which, respectively, can be interpreted as:

$$\begin{array}{}\text{(15)}& \mathrm{\Delta }v\left(\mathrm{\Delta }x\right)\stackrel{d}{=}\mathit{\epsilon }(\mathrm{\Delta }x{)}^{a_{\mathrm{h}}}\mathrm{\Delta }{x}^{H_{\mathrm{h}}};\mathrm{\Delta }v(\mathrm{\Delta }z\left)\stackrel{d}{=}\mathit{\varphi }\right(\mathrm{\Delta }z{)}^{a_{\mathrm{v}}}\mathrm{\Delta }{z}^{H_{\mathrm{v}}},\end{array}$$

with scaling exponents $a_{\mathrm{h}}=\mathrm{1}/\mathrm{3}$ for the kinetic energy flux density ε and $a_{\mathrm{v}}=\mathrm{1}/\mathrm{5}$ for the buoyancy force flux density ϕ. Due to the fact that eddies can be defined as structures having a given velocity shear Δv, we have, in fact, the following balance between both flux densities:

$$\begin{array}{}\text{(16)}& \mathit{\epsilon }\left(\mathrm{\Delta }x{)}^{a_{\mathrm{h}}}\mathrm{\Delta }{x}^{H_{\mathrm{h}}}\stackrel{d}{=}\mathit{\varphi }\right(\mathrm{\Delta }z{)}^{a_{\mathrm{v}}}\mathrm{\Delta }{z}^{H_{\mathrm{v}}}\end{array}$$

and finally:

$$\begin{array}{}\text{(17)}& \mathit{\epsilon }\left(\mathrm{\Delta }x{)}^{a_{\mathrm{h}}}\stackrel{d}{=}\mathit{\varphi }\right(\mathrm{\Delta }z{)}^{a_{\mathrm{v}}},\mathrm{with}:\mathrm{\Delta }{x}^{H_{\mathrm{h}}}\approx \mathrm{\Delta }{z}^{H_{\mathrm{v}}},\end{array}$$

due to the respective scaling of the horizontal and vertical eddy sizes defined by the scale change T_λ (Eq. 12).

This confirms that ϕ and ε are both sides of the same coin (Eq. 11). In summary, what is only needed is an anisotropic cascade defined by an anisotropic scale compatible (e.g. Eq. 14) with the anisotropic scale change T_λ (Eq. 12), the intermittency being introduced either by the flux density ε or ϕ (Eq. 15), which respect their equivalence (Eq. 17).

Practical difficulties only occur with discrete scales because they cannot easily deal with arbitrary H_z.

Discrete scales and an ad-hoc approximation

To get around these difficulties, it was tentatively proposed to consider the following 3D model:

$$\begin{array}{}\text{(18)}& \begin{array}{l}{v}_x(x+\mathrm{\Delta }x,y+\mathrm{\Delta }y,z+\mathrm{\Delta }z,t)\\ =v_x(x,y,z,t)+c_{{\mathit{ϵ}}_x}{\mathit{ϵ}}_x(x+\mathrm{\Delta }x,y+\mathrm{\Delta }y,t{)}^{a_{\mathrm{h}}}\mathrm{\Delta }{l}^{H_{\mathrm{h}}}\\ +c_{{\mathit{\theta }}_x}{\mathit{\theta }}_x(z+\mathrm{\Delta }z,t{)}^{a_{\mathrm{v}}}\mathrm{\Delta }{z}^{H_{\mathrm{v}}}\\ v_y(x+\mathrm{\Delta }x,y+\mathrm{\Delta }y,z+\mathrm{\Delta }z,t)\\ =v_y(x,y,z,t)+c_{{\mathit{ϵ}}_y}{\mathit{ϵ}}_y(x+\mathrm{\Delta }x,y+\mathrm{\Delta }y,t{)}^{a_{\mathrm{h}}}\mathrm{\Delta }{l}^{H_{\mathrm{h}}}\\ +c_{{\mathit{\theta }}_y}{\mathit{\theta }}_y(z+\mathrm{\Delta }z,t{)}^{a_{\mathrm{v}}}\mathrm{\Delta }{z}^{H_{\mathrm{v}}}\\ v_z(x+\mathrm{\Delta }x,y+\mathrm{\Delta }y,z+\mathrm{\Delta }z,t)\\ =v_z(x,y,z,t)+c_{{\mathit{ϵ}}_z}{\mathit{ϵ}}_z(x+\mathrm{\Delta }x,y+\mathrm{\Delta }y,t{)}^{a_{\mathrm{h}}}\mathrm{\Delta }{l}^{H_{\mathrm{h}}}\\ +c_{{\mathit{\theta }}_z}{\mathit{\theta }}_z(z+\mathrm{\Delta }z,t{)}^{a_{\mathrm{v}}}\mathrm{\Delta }{z}^{H_{\mathrm{v}}},\end{array}\end{array}$$

where $\mathrm{\Delta }l=\sqrt{\mathrm{\Delta }{x}^{\mathrm{2}}+\mathrm{\Delta }{y}^{\mathrm{2}}}$ is the (classical), horizontal distance from the anenometer and θ is the buoyancy force flux but is supposed independent of ε contrary to ϕ (Eq. 17). Unfortunately, this model suffers from a number of basic problems:

• Equalities in distribution are replaced by deterministic equalities, which oversimplify and trivialise the dynamics for each realisation.

• The flux density ε is defined only along three directions instead along all directions (Eq. 11), whose corresponding components (ε_x, ε_y, ε_z) are supposed to be independent, which is not a tenable assumption.

• The introduction of θ independent of ε is purely ad hoc to additively(!) introduce a second isotropic scaling (!), and therefore the change of scaling is reduced to a linear cross-over.

• The scale of the flux densities ε and θ is implicitly taken as the scale of the simulation resolution instead of the pair separation scale (Eq. 15). The resulting mismatch between both scales will introduce a statistical bias.

• Moreover, the density values are arbitrarily taken at the locations ($x+\mathrm{\Delta }x,y+\mathrm{\Delta }y$ for ε) and z+Δz for θ. This introduces a dissymmetry that is not without consequences.

• All time steps are fully independent, with the exception of the empirical velocity values $v(x,y,z,t)$ measured by the 3D sonic anemometer.

The limitations of this model are, therefore, extremely strong. Many of them would have been resolved with the help of the scalar anisotropic cascades recalled above (see the previous paragraph). But to fully overcome them would require us to consider their extension to vector fields (Schertzer and Tchiguirinskaia, 2015, 2020). This remains outside the scope of this paper, and this simplistic approach is used for a sort of exploratory step focused on the interaction of drops and a given velocity field. The results obtained, therefore, have to be critically examined, since many of them are sensitive to the aforementioned oversimplifications.

4.1.3 Practical implementation

The fields ϵ (i.e. the ones for the horizontal shift) are simulated in space time with a size $\mathrm{729}\times \mathrm{729}\times \mathrm{64}$ using discrete UM cascades and Eq. (10) to obtain either positive or negative values. The number 729 refer to space, while 64 refers to time. Before stating the physical resolution of the field, it should be explained that a simple anisotropy between space and time is accounted through a scaling anisotropy coefficient H_t. Within such framework, when the spatial scale of the data is changed by a ratio of λ_xy, the temporal scale should be changed by a factor of ${\mathit{\lambda }}_t={\mathit{\lambda }}_{xy}^{H_t}$; H_t is expected to be equal to $\mathrm{1}/\mathrm{3}$ (Marsan et al., 1996), hence when the spatial scale is multiplied by 3, the temporal scale should be multiplied by 2 (i.e. ${\mathrm{3}}^{\mathrm{1}-\mathrm{1}/\mathrm{3}}\approx \mathrm{2.08}$) (Biaou et al., 2005; Gires et al., 2014). These fields are assumed to cover an area of the size 40 km × 40 km × 1024 s, which is needed for drift of 0.1 mm drops during their fall when wind is strong. This means that a voxel is of size 53 m × 53 m × 16 s.

The fields θ (i.e. those for the vertical shift) are of size 512×64 covering a physical are of 1600 m × 1024 s, meaning a pixel size 3 m × 16 s. All UM fields are simulated with α=1.7, C₁=0.2, as in the previous section.

Finally, at any point $(x,y,z,t)$ a bi or tri-linear interpolation is implemented to obtain the value of the field from the nearest points. The value of the prefactor were set to $c_{{\mathit{ϵ}}_x}=c_{{\mathit{ϵ}}_y}=\mathrm{0.3}$, $c_{{\mathit{ϵ}}_z}=\mathrm{0.1}$ and $c_{{\mathit{\theta }}_x}=c_{{\mathit{\theta }}_y}=c_{{\mathit{\theta }}_z}=\mathrm{0.01}$ through an heuristic approach of trial and error to get some realistic fluctuations. In the future, it would obviously be necessary to tune them to local wind properties. However, such tuning is outside the scope of this section, which aims more at being a proof of concept.

4.2 Illustration

In order to illustrate the suggested process, let us consider a 0.5 mm drop during the low wind event. Its initial position is $(\mathrm{0},\mathrm{0},\mathrm{1500})$ in m. It is “dropped” with no horizontal velocity and a vertical velocity equal to that of its terminal. The anemometer is assumed to be located at $(\mathrm{0},\mathrm{0},\mathrm{100})$ m. Then Eq. (5) is implemented. At each time step, the local wind is assessed using the methodology described in the previous paragraphs.

The actual total wind perceived by the drop (i.e. the input in Eq. 5) is recorded and displayed in last row of Fig. 7. It corresponds to the sum of the wind from the anemometer (the first row in Fig. 7) plus a wind shift field (the middle row in Fig. 7). This leads to a given trajectory in space, which is shown in Fig. 8. The projected trajectory on the planes (x,y) and (x,z) is also shown. This trajectory exhibits a nonlinear complex pattern, which results from the turbulent nature of the wind. It should be mentioned that not accounting for oblateness affects not only the vertical velocity but the whole trajectory as well. Depending on the model (spheres or oblate drops), shifts of more than few hundred metres are found even for large drops for some events.

Figure 7 (Top) Temporal evolution with 0.01 s time steps of the wind data from 3D sonic anemometer, (middle) the wind shift and (bottom) the total wind perceived by the 0.5 mm drop falling from the position $(\mathrm{0},\mathrm{0},\mathrm{1500})$ during the low wind event. Bottom panel is actually the wind input used to obtain the trajectory of Fig. 8.

Figure 8 Trajectory (solid line) of a 0.5 mm drop in a turbulent wind field for the low wind event. The dotted lines correspond to the trajectory projected on the (x,z) and (y,z) plane.

4.3 Sensitivity to the wind shift field

The process to generate an estimation of a 3D wind field is actually stochastic through the UM fields ϵ and θ used in Eq. (18). In this section, the sensitivity to the given realization of the process is discussed. In order to achieve that, 10 wind samples are generated with the same input parameters (described in Sect. 4.1) and the corresponding trajectories for drops of size 0.5, 1, 2 and 3 mm are computed.

The projected trajectories for the low wind event, are displayed in Fig. 9a and b. The position of the dropa when they reach the ground is shown in Fig. 9c. The spread of the drops strongly depends on their size, with a decrease as the drop size increases. Indeed, Δx (x_max−x_min) is equal to 1238, 591, 415 and 404 m for drops of sizes 0.5, 1, 2 and 3 mm, respectively. For Δy, the values are 2123, 897, 461 and 322 m. Such a decrease is due to a combination of the fact that smaller drops are more subject to wind fluctuations (Sect. 3) and that they spend more time in the atmosphere (Sect. 2) before they reach the ground. Similar trends are retrieved for the strong wind event (Fig. 10) with a stronger absolute shift. In that case, the values for Δx are 890, 868, 500 and 448 m, respectively. For Δy, they are 2382, 1003, 568 and 448 m.

Figure 9 (a, b) Trajectories of drops of various sizes falling from the $(\mathrm{0},\mathrm{0},\mathrm{1500})$ position, projected on the (x,z) and (y,z) plane, respectively. The low wind event is used. The ten curves correspond to different realizations of the wind shift field. (c) Position of the ground impact of the various drops.

Figure 10 The same as Fig. 9 but for the strong wind event.

4.4 Illustration of impact on rainfall retrieval with weather radars

In this last section, initial investigations toward understanding the consequence of previous work on quantitative rainfall measurement with weather radars are carried out. Indeed, weather radar measure rainfall at a given altitude, while hydro-meteorologists are interested in rainfall at ground level. During their fall, significant shifts can occur. In order to study it, the following process is implemented. For 5 min, one rainfall drop is dropped every 15 s from a random position within a voxel of size 100 m centred on $(\mathrm{0},\mathrm{0},\mathrm{1450})$ m. Hence it covers a total duration of 5 min. The trajectories and positions on the ground of the drops is then studied. This enables us to basically mimic the measurement of a weather radar at its typical gate size and temporal resolution.

Figure 11 displays the trajectories and ground impact location in the case of the low wind event for a given realisation of the stochastic wind shift. A shift of more than 1 km for small drops is found, and more than 300 m for 3 mm drops. As noted in the previous section, the spread of drops at ground level tends to decrease with increasing drop size. Indeed, Δx (x_max−x_min) is equal to 206, 119, 165 and 121 m for drop of sizes 0.5, 1, 2 and 3 mm, respectively. For Δy, the values are 184, 134, 112 and 88 m. For the strong wind event (Fig. 12), shifts of more than 6 and 1.5 km are reported for drops of size 0.5 and 3 mm respectively. Similar results as for the low wind event are found with regards to the spread. The corresponding figures are of 665, 450, 284 and 307 m for Δx, and of 819, 664, 429 and 420 m for Δy.

Figure 11 (a, b) Trajectories of drops of various sizes falling from a 100 m cubic voxel centred on a $(\mathrm{0},\mathrm{0},\mathrm{1450})$ position, projected on the (x,z) and (y,z) plane, respectively. The low wind event is used. For each size, the 20 drops are dropped every 15 s (hence over a total duration of 5 min). A single realization of the wind shift field is used. (c) Position of the ground impact of the various drops.

Figure 12 The same as in Fig. 11 but for the strong wind event.

As was previously pointed out, this spread is due to the fact that smaller drops spend more time in the atmosphere and are more sensitive to wind fluctuations. Indeed, the duration of a fall from 1500 m to the ground at 0 m is equal to 716, 378, 238 and 192 s for drops of size 0.5, 1, 2 and 3 mm, respectively. Given that high-resolution radar pixels are typically of the size of a few hundred metres, one should note that drops within a given voxel at measurement height can reach the ground within an area of size 3 km × 6 km. Even within a drop diameter class, shifts are of few radar pixels. Given that drop size distribution also varies, such a shift can significantly affect rainfall retrieval, even in low wind conditions.

5 Conclusions

In this paper, we have aimed for a better understanding of the behaviour of individual rainfall drops falling from typically 1500 m. In a first step, we developed a new approach to compute the drag coefficient accounting for drop oblateness and findings in fluid mechanics. This was validated for drops of equivolumic size of up to 4 mm through the comparison between the retrieved terminal fall velocity and the commonly used formula.

Then the temporal evolution of the horizontal drop velocity under turbulent wind constraints was studied. It appears that multifractal features of the input wind are simply transferred to drop velocity with an additional fractional integration and slight time shift. The UM parameters α and C₁ are basically conserved, while H is increased. The increase ranges from 0.1 for 0.1 mm size drops to 0.8 for 1–1.5 mm size drops. It remains rather constant for larger drops.

Finally, the trajectories of drops of various sizes falling form 1500 m was studied as a proof of concept. For this, 100 Hz anemometer data was used, and an approach to simulate realistic fluctuations of wind in space was developed. It notably enables to analyse how drop shift during their fall between their location measurement by weather radars and ground impact. For a strong wind event, drops located within a radar gate in altitude for 5 min are spread on the ground over an area of a few kilometres. The spread for drops of a given diameter is found to cover a few radar pixels.

In order to further explore the consequences of these findings on quantitative rainfall estimation with weather radars, further investigations are needed. More precisely, (i) the model to simulate wind fluctuations should be improved, notably to use vector simulations and tune the prefactors according to local wind conditions; (ii) space-time outputs of numerical weather prediction models could also be tested to retrieve wind fields; (iii) the actual drop size distribution should be used to better assess the impact for the ground estimation of precipitation, which implies making some simulations for a much larger number of drops; (v) a longer period of time should be tested to investigate where the water volume (i.e. all the drops) of a given radar gate fall during an event. For the two last points, data is available within the RW-Turb project. Such step would then need to be repeated over various radar gates to derive updated radar maps. Given the limited computation power that will not allow us to simulate the trajectories of all the drops, some statistical behaviour according to each radar gate and wind conditions would need to be designed and then computed. It should also be stressed that only individual drops are currently being handled. This means that the methodology developed does not account for either collision, aggregation between drops or for breakup. Such processes are also known to affect drop velocities by changing their size and shape. Future investigations should also aim at accounting for them. Finally, it should also be stressed that the method developed stochastically simulates wind fluctuations at small scales. This means that the output will not be a deterministic radar measurement but an ensemble of possible realistic outputs, out of which a probability distribution could be derived. Such a probabilistic approach is discussed in Kirstetter et al. (2015) with a focus on intrinsic radar uncertainties and not wind drift.