Articles | Volume 21, issue 17
Atmos. Chem. Phys., 21, 13119–13130, 2021
https://doi.org/10.5194/acp-21-13119-2021
Atmos. Chem. Phys., 21, 13119–13130, 2021
https://doi.org/10.5194/acp-21-13119-2021

Research article 03 Sep 2021

Research article | 03 Sep 2021

Statistical properties of a stochastic model of eddy hopping

Statistical properties of a stochastic model of eddy hopping
Izumi Saito, Takeshi Watanabe, and Toshiyuki Gotoh Izumi Saito et al.
  • Department of Physical Science and Engineering, Nagoya Institute of Technology, Gokiso-cho, Showa-ku, Nagoya 466-8555, Japan

Correspondence: Izumi Saito (izumi@gfd-dennou.org)

Abstract

Statistical properties are investigated for the stochastic model of eddy hopping, which is a novel cloud microphysical model that accounts for the effect of the supersaturation fluctuation at unresolved scales on the growth of cloud droplets and on spectral broadening. Two versions of the model, the original version by Grabowski and Abade (2017) and the second version by Abade et al. (2018), are considered and validated against the reference data taken from direct numerical simulations and large-eddy simulations (LESs). It is shown that the original version fails to reproduce a proper scaling for a certain range of parameters, resulting in a deviation of the model prediction from the reference data, while the second version successfully reproduces the proper scaling. In addition, a possible simplification of the model is discussed, which reduces the number of model variables while keeping the statistical properties almost unchanged in the typical parameter range for the model implementation in the LES Lagrangian cloud model.

1 Introduction

The purpose of the present paper is to investigate the statistical properties of the stochastic model of eddy hopping proposed by Grabowski and Abade (2017). This stochastic model, referred to hereinafter as the eddy-hopping model, was developed in order to account for the effect of the supersaturation fluctuation at unresolved (subgrid) scales on the growth of cloud droplets by the condensation process. In a turbulent cloud, cloud droplets arriving at a given location follow different trajectories and thus experience different growth histories, which leads to significant spectral broadening. This mechanism, referred to as the stochastic condensation theory, has been investigated since the early 1960s by a number of researchers (mostly Russian; see Sedunov1974; Clark and Hall1979; Korolev and Mazin2003), but the importance of this mechanism was later reinforced by Cooper (1989) and Lasher-Trapp et al. (2005). For this mechanism, Grabowski and Wang (2013) emphasized the importance of large-scale eddies (turbulent eddies with scales not much smaller than the cloud itself) and proposed the concept of large-eddy hopping. Grabowski and Abade (2017) formulated this concept and developed the eddy-hopping model. Abade et al. (2018) extended the model by introducing a term accounting for the relaxation of supersaturation fluctuations due to turbulent mixing. For clarity, we hereinafter refer to the model by Grabowski and Abade (2017) as the original version, and the model by Abade et al. (2018) as the second version. For the following study using the eddy-hopping model, readers are referred to Thomas et al. (2020).

It should be noted that the turbulent entrainment mixing is another important mechanism for the supersaturation fluctuation generation other than the stochastic condensation and that the effects of the turbulent entrainment mixing are not included in the eddy-hopping model considered in the present study. Abade et al. (2018) investigated the effects of the turbulent entrainment mixing and entrained CCN activation by using the entraining parcel model.

In the present paper, we take a rather theoretical approach to obtain various statistical properties of the eddy-hopping model, such as the variance, covariance, and auto-correlation function of the supersaturation fluctuation. These statistical properties are used to validate the model against the reference data taken from direct numerical simulations (DNSs) and large-eddy simulations (LESs). We show that the original version of the eddy-hopping model fails to reproduce a proper scaling for a certain range of parameters, resulting in the deviation of the model prediction from the reference data, while the second version successfully reproduces the proper scaling. We show how the relaxation term introduced by Abade et al. (2018) leads to this improvement. We also discuss the possibility of simplification of the model, which reduces the number of model variables while keeping the statistical properties almost unchanged in the typical parameter range for the model implementation in the LES Lagrangian cloud model.

The remainder of the present paper is organized as follows. Section 2 describes the governing equations of the original version. Section 3 presents a theoretical analysis and numerical experiments and demonstrates the improper scaling in the model prediction by the original version. Section 4 describes the second version. Finally, Sect. 5 discusses the possibility of simplification of the model.

2 Governing equations

The original version of the eddy-hopping model proposed by Grabowski and Abade (2017) consists of the following two evolution equations. First, the fluctuation of the vertical velocity of turbulent flow at the droplet position, w(t), is modeled by the Ornstein–Uhlenbeck process:

(1) w ( t + δ t ) = w ( t ) e - δ t / τ + 1 - e - 2 δ t τ σ w ψ ,

where δt is the time increment, ψ is a Gaussian random number with zero mean and unit variance drawn every time step, σw is the standard deviation of w, and τ is the integral time, or the large-eddy turnover time of the turbulent flow. Here, σw and τ are used as the model parameters. Second, the supersaturation fluctuation at the droplet position, S(t), is governed by

(2) d S d t = a 1 w - S τ relax .

Here, the first term on the right-hand side represents the effect of adiabatic cooling and warming due to air parcel ascent and descent caused by the vertical velocity w(t). The parameter a1 has the unit of a scalar gradient. The second term on the right-hand side represents the effect of condensation or evaporation of droplets. The timescale τrelax is referred to as the phase relaxation time and is inversely proportional to the average of the number density and radius of the droplets (Politovich and Cooper1988; Korolev and Mazin2003; Kostinski2009; Devenish et al.2012).

Equation (1) can also be written as the following derivative form (Pope2000):

(3) d w d t = - 1 τ w ( t ) + F w ( t ) .

Here, the term Fw(t) is statistically independent of S and obeys the Gaussian random process that has zero mean and two-time covariance defined by

(4) F w ( t ) F w ( s ) = 2 σ w 2 τ δ ( t - s ) ,

where the angle brackets indicate an ensemble average and δ( ) is the Dirac delta function. In the following theoretical analysis, Eqs. (3) and (2) are used as the governing equations of the original version.

3 Statistical properties of the original version

We now obtain the analytical expression for the standard deviation of the supersaturation fluctuation, σS, in a statistically steady state. Starting from Eqs. (3) and (2), the result is provided in Eq. (13).

First, multiplying Eq. (3) by S and taking an ensemble average, we obtain

(5) S d w d t = - 1 τ w S

because of statistical independence SFw=0. Second, multiplying Eq. (2) by w and taking an ensemble average, we obtain

(6) w d S d t = a 1 w 2 - 1 τ relax w S .

Summing Eqs. (5) and (6), we obtain

(7) d d t w S = a 1 w 2 - 1 τ relax w S - 1 τ w S .

Next, we consider a statistically steady state. Since an ensemble-averaged variable does not change in time (d/dt=0) and w2=σw2 in the statistically steady state, we obtain the flux of the supersaturation in the vertical direction as follows:

wS=a11τ+1τrelax-1σw2(8)=a1(1+Da)-1τσw2,

where Da is the Damköhler number (Shaw2003) defined as

(9) D a = τ τ relax .

Next, multiplying Eq. (2) by S and taking an ensemble average, we obtain

(10) 1 2 d d t S 2 = a 1 w S - 1 τ relax S 2 .

In the statistically steady state, we have

(11) σ S 2 = S 2 = a 1 τ relax w S .

Combining Eqs. (8) and (11), we obtain

σS2=a1τrelaxa1(1+Da)-1τσw2(12)=a12(1+Da)-1τrelaxτσw2,

or equivalently,

(13) σ S = ( 1 + D a ) - 1 / 2 D a - 1 / 2 a 1 τ σ w .

Here, σS in Eq. (13) has two important asymptotic forms, as shown below:

  1. Large scale limit.
    For τ→∞ (or equivalently, Da→∞, L→∞, where L=σwτ is the integral scale), σS in Eq. (13) is approximated as

    σSa1Da-1/2τrelax1/2τ1/2σw(14)=a1τrelaxσw.

    For the case of a constant dissipation rate of turbulent kinetic energy ε, σwL1/3 (see Appendix B), and we have the following scaling:

    (15) σ S L 1 / 3 .
  2. Small scale limit.
    For τ→0 (or equivalently, Da→0, L→0), σS in Eq. (13) is approximated as

    (16) σ S a 1 τ relax 1 / 2 τ 1 / 2 σ w .

    For the case of a constant dissipation rate of turbulent kinetic energy ε, σwL1/3 and τL2/3 (see Appendix B), and we have the following scaling:

    (17) σ S L 2 / 3 .

The above asymptotic forms of σS in the two limits can be validated through comparison with the result of the scaling argument by Lanotte et al. (2009). From their argument, we should have σSa1τrelaxσw for the large scale limit, which is consistent with Eq. (14). On the other hand, we should have σSa1τσw for the small scale limit, which is inconsistent with Eq. (16). Therefore, the original version given by Eqs. (3) and (2) does not reproduce the proper scaling for the small scale limit.

https://acp.copernicus.org/articles/21/13119/2021/acp-21-13119-2021-f01

Figure 1 Standard deviation of the supersaturation fluctuation σS in the statistically steady state obtained from the analytical expression given by Eq. (13) (orange curve) and the results of our numerical integration of the original version of the eddy-hopping model (blue squares). The horizontal axis is the integral length L. The black dots indicate the reference data taken from direct numerical simulations and large-eddy simulations by Thomas et al. (2020). The red triangles indicate the results of the numerical integration of the original version reported by Thomas et al. (2020). The range of L and σS for the panel is the same as in Fig. 10 in Thomas et al. (2020). The three short black lines indicate slopes of 1, 2/3, and 1/3.

Download

Figure 1 compares the scale dependence of σS for the analytical expression given by Eq. (13) (orange curve) with the results of the numerical integration of the original version given by Eqs. (1) and (2) (blue squares). Here, numerical integration is conducted in the same manner as that by Thomas et al. (2020) (Sect. 5 of their study), except that the integration time is increased from 6τ to 10τ (see Appendix A for details). After the integration time of 10τ, all of the experimental results achieved a statistically steady state and agreed with the theoretical curve (compare the orange curve and the blue squares). As expected based on the analysis, the theoretical curve shows the scaling σSL1/3 for large scales (approximately L>101m) and the improper scaling σSL2/3 for small scales (approximately L<10-1m). These results are contrary to the results of DNSs and LESs (scaled-up DNSs) conducted by Thomas et al. (2020) (black dots in Fig. 1), which show the proper scalings both for large and small scales (σSL1/3 and L1, respectively).

Note that Fig. 1 also shows the results of the numerical integration of the original version reported by Thomas et al. (2020) (red triangles), and their results disagree with the results of the present study. A possible reason for this discrepancy might be that their results did not achieve a statistically steady state. For details, see Appendix C.

The original version of the eddy-hopping model given by Eqs. (1) and (2) shows the improper scaling for small scales because of the assumption made in the formulation of the model. Originally, Eq. (2) (corresponding to Eq. 8 in Grabowski and Abade2017) was formulated under the assumption of large scales (or Da≫1), since this assumption usually holds for typical situations in atmospheric clouds. Thus, it is reasonable that the original version does not reproduce the proper scaling for small scales.

4 Statistical properties of the second version

We next consider the second version of the eddy-hopping model by Abade et al. (2018), which is written as follows:

(18)w(t+δt)=w(t)e-δt/(c1τ)+1-e-2δt(c1τ)σwψ,(19)dSdt=a1w-S(c2τrelax)-S(c1τ).

Note that, for subsequent use, we write the governing equations in a slightly generalized form by introducing two parameters c1 and c2. The second version by Abade et al. (2018) has c1=c2=1.

The important change introduced by Abade et al. (2018) into the original version is the term proportional to -S/τ in Eq. (19). Physically, this term represents the damping effect on S due to turbulent mixing (eddy diffusivity). This type of term is commonly included in stochastic models used in cloud turbulence research (Sardina et al.2015, 2018; Chandrakar et al.2016; Siewert et al.2017; Saito et al.2019a). The timescale of the damping effect due to turbulent mixing is characterized by the integral time τ, whereas that due to condensation or evaporation of cloud droplets is characterized by the phase relaxation time τrelax. The relative importance of these two effects is characterized by the Damköhler number (Da=τ/τrelax), where the damping effect due to turbulent mixing is dominant for Da≪1 (corresponding to small scales). Below we show that the term -S/τ plays an essential role in reproducing the proper scaling.

Note that it is possible to further extend the second version by additionally introducing the Wiener process term representing small-scale fluctuations or mixing into Eq. (19) for the supersaturation fluctuation. Readers are referred to Paoli and Shariff (2009) and Sardina et al. (2015) for the Langevin model including such terms, and also to Chandrakar et al. (2021), who implemented such a Langevin model into the LES Lagrangian cloud model and investigated the impact of entrainment mixing and turbulent fluctuations on droplet size distributions in a cumulus cloud. In the present study, however, we focus on statistical properties of the second version with Eqs. (18) and (19).

Applying the analytical procedure described in Sect. 3 to the second version given by Eqs. (18) and (19), we first obtain

wS=a12c1τ+1c2τrelax-1σw2(20)=c1a1[2+(c1/c2)Da]-1τσw2,

instead of Eq. (8). Next, instead of Eq. (11), we have

σS2=S2=a11c1τ+1c2τrelax-1wS(21)=c1a1[1+(c1/c2)Da]-1τwS.

Finally, the analytical expression corresponding to Eq. (13) is

(22) σ S = [ 1 + ( c 1 / c 2 ) D a ] - 1 / 2 [ 2 + ( c 1 / c 2 ) D a ] - 1 / 2 c 1 a 1 τ σ w .

Asymptotic forms of σS in Eq. (22) for the large and small scale limits are, respectively, given as follows:

  1. Large scale limit.
    For τ→∞ (or equivalently, Da→∞, L→∞), σS in Eq. (22) is approximated as

    σSc2a1Da-1τσw(23)=c2a1τrelaxσw.

    For the case of a constant dissipation rate of turbulent kinetic energy ε, we have

    (24) σ S L 1 / 3 .
  2. Small scale limit.
    For τ→0 (or equivalently, Da→0, L→0), σS in Eq. (22) is approximated as

    (25) σ S 2 - 1 / 2 c 1 a 1 τ σ w ,

    which indicates that τrelax1/2 in Eq. (16) has been replaced by τ1/2 by introducing the term -S/τ in Eq. (19). For the case of a constant dissipation rate of turbulent kinetic energy ε, we have

    (26) σ S L .

Therefore, the second version successfully reproduces asymptotic forms σSa1τrelaxσw and a1τσw for the large and small scale limits, respectively, which are both consistent with the result of the scaling argument by Lanotte et al. (2009).

https://acp.copernicus.org/articles/21/13119/2021/acp-21-13119-2021-f02

Figure 2 Standard deviation of the supersaturation fluctuation σS in the statistically steady state obtained from the analytical expression given by Eq. (22) for the second version (c1=c2=1, orange curve) and the results of our numerical integration using the second version given by Eqs. (18) and (19) (c1=c2=1, green diamonds). The dashed line indicates the analytical expression given by Eq. (13) for the original version. The two short black lines indicate slopes of 1 and 1/3. The black dots and the axes of the panel are the same as in Fig. 1.

Download

Figure 2 (orange curve) shows the theoretical curve given by Eq. (22) for the second version (c1=c2=1). The second version reproduces the proper scalings both for large and small scales (σSL1/3 and L1, respectively) and demonstrates better agreement with the reference data (black dots) than the original version for L<100 m. Figure 2 (green diamonds) also shows the results of the numerical integration of the second version, which agree with the theoretical curve (orange curve), as expected. Here, the numerical integration was conducted in the same manner as in the previous section (see Appendix A for details).

https://acp.copernicus.org/articles/21/13119/2021/acp-21-13119-2021-f03

Figure 3 Standard deviation of the supersaturation fluctuation σS in the statistically steady state from the analytical expression (22) for the second version. The dashed curve is for c1=c2=1, and the solid curve is for c1=0.746 and c2=1.28. The black dots and the axes of the panel are the same as in Fig. 1.

Download

Although improved, the second version still slightly over- and underestimates the supersaturation fluctuations for L<3×10-1 and L>2×100 m, respectively, as shown in Fig. 2 (compare the orange curve with black dots). This deviation from the reference data can be further reduced by adjusting two parameters c1 and c2 in Eqs. (18) and (19). The analytical expression (22) and its asymptotic forms (23) and (25) show how c1 and c2 work. These types of parameters are not new. For example, a parameter corresponding to c1 is commonly used in the Langevin stochastic equation in turbulence research (Sawford1991; Marcq and Naert1998). Formally, the inverse of c1 is referred to as the drift coefficient, and the coefficients for the velocity and scalar equations should be distinguished. However, we treat these coefficients as the same parameter in Eqs. (18) and (19) for simplicity. On the other hand, the importance of a parameter corresponding to c2 has been demonstrated in a recent study on turbulence modulation by particles (Saito et al.2019b).

Here, we do not consider any physical meaning for c1 and c2 and use them just as tuning parameters. Two parameters c1 and c2 are determined by comparing the theoretical curve given by Eq. (22) with the reference data taken from DNSs and LESs in Thomas et al. (2020). The best fit is given by c1=0.746 and c2=1.28. Figure 3 (solid curve) shows the theoretical curve given by Eq. (22) with these values of c1 and c2, which agrees almost perfectly with the reference data (black dots). Although the improvement from the second version with c1=c2=1 is slight, this result shows that the eddy-hopping model can be easily tuned to reproduce the reference data almost perfectly.

5 Possibility of simplification of the model

Finally, we discuss the possibility of simplification of the eddy-hopping model. Here, our discussion is based on the second version given by Eqs. (18) and (19), but the same argument also applies to the original version given by Eqs. (1) and (2).

https://acp.copernicus.org/articles/21/13119/2021/acp-21-13119-2021-f04

Figure 4 Auto-correlation functions in the statistically steady state for the simplified model (B(t) in Eq. 35: blue dashed curve) and the second version (A(t) in Eq. 28: red solid curve). Here, c1=c2=1. The parameters for each panel are as follows: (a) L=10-2 m, τ=0.447 s, τ0=0.844 s, Da=0.127; (b) L=10-1 m, τ=2.08 s, τ0=3.38 s, Da=0.591; (c) L=100 m, τ=9.63 s, τ0=12.2 s, Da=2.74; (d) L=101 m, τ=44.7 s, τ0=48.0 s, Da=12.7; and (e) L=102 m, τ=208 s, τ0=211 s, Da=59.1. The phase relaxation time is fixed to τrelax=3.513 s. The horizontal axis is the time t normalized by the auto-correlation time τ0 for each case. The parameter τ is determined from the integral length L based on the setting for the numerical experiment described in Appendices A and B.

Download

The eddy-hopping model consists of two evolution equations for the supersaturation and vertical velocity fluctuations, S and w respectively, and these two variables fluctuate randomly according to the Ornstein–Uhlenbeck process. However, if we have S that fluctuates with a proper amplitude and auto-correlation function, then we do not need the evolution equation for w, because only S is used in the growth equation of the droplet size. As described in Sect. 4, we obtained an analytical expression for σS, i.e., the standard deviation of the supersaturation fluctuation in the statistically steady state given by Eq. (22). On the other hand, the auto-correlation function for S in a statistically steady state can also be obtained analytically. The derivation is described in Appendix D. The result is given in Eq. (D14) and is as follows:

(27)A(t)=S(t+t0)S(t0)S(t0)S(t0)(28)=τ1τ1-τ2e-t/τ1-τ2τ1-τ2e-t/τ2,

where τ1 and τ2 are, respectively, defined as

(29) τ 1 = c 1 τ , and τ 2 = 1 c 1 τ + 1 c 2 τ relax - 1 = c 1 τ 1 + c 1 c 2 D a - 1 .

We can also obtain the auto-correlation time for S by time integration of A(t) (see Eq. D16 in Appendix D), which is given as

(30) τ 0 = τ 1 + τ 2 .

The auto-correlation function A(t) in Eq. (28) and the auto-correlation time τ0 in Eq. (30), with τ1 and τ2 defined by (29), have important asymptotic forms in two limits. First, for the large scale limit (Da→∞), the asymptotic forms of A(t) and τ0 are given by

(31) lim D a A ( t ) = e - t / ( c 1 τ ) and lim D a τ 0 = c 1 τ ,

respectively. Second, for the small scale limit (Da→0), the asymptotic forms of A(t) and τ0 are given by

(32) lim D a 0 A ( t ) = 1 + t ( c 1 τ ) e - t / ( c 1 τ ) and lim D a 0 τ 0 = 2 c 1 τ ,

respectively.

https://acp.copernicus.org/articles/21/13119/2021/acp-21-13119-2021-f05

Figure 5 Time evolutions of the supersaturation fluctuations obtained from the numerical integration of the original version (Eqs. 1 and 2, dashed red curve), the second version (Eqs. 18 and 19, dotted blue curve), and the simplified model (Eq. 33, solid green curve). The parameters for each panel are as follows: (a) L=10-2 m, τ=0.447 s, Da=0.127; (b) L=100 m, τ=9.63 s, Da=2.74; (c) L=102 m, τ=208 s, Da=59.1. The phase relaxation time is fixed to τrelax=3.513 s. Results are shown for the time range 0(t/τ-10)10 (or equivalently, 10τt20τ), where all cases are already in a statistically steady state. The numerical integration was conducted in the same manner as described in Appendix A. All results were obtained by using the same random number series.

Download

Based on analytical expressions for the fluctuation amplitude and the auto-correlation function for S (Eqs. 22 and 30, respectively), a simplified version of the eddy-hopping model is defined as follows:

(33) S ( t + δ t ) = S ( t ) e - δ t / τ 0 + 1 - e - 2 δ t τ 0 σ S ψ ,

where σS and τ0 are given by Eqs. (22) and (30), respectively. Note that the simplified model given by Eq. (33) is a single-equation model, as compared to the two-equation model given by Eqs. (18) and (19) before the simplification. The auto-correlation function for S in the simplified model given by Eq. (33) is as follows:

(34)B(t)=S(t+t0)S(t0)S(t0)S(t0)(35)=e-t/τ0,

which has the following two asymptotic forms. First, for the large scale limit (Da→∞),

(36) lim D a B ( t ) = e - t / ( c 1 τ ) ,

which agrees with the corresponding asymptotic form given by Eq. (31) for the second version. Second, for the small scale limit (Da→0),

(37) lim D a 0 B ( t ) = e - t / ( 2 c 1 τ ) ,

which disagrees with the corresponding asymptotic form given by Eq. (32) for the second version.

Figures 4a through e compare the auto-correlation function for the simplified model (B(t) in Eq. 35: blue dashed curve) and that for the second version (A(t) in Eq. 28: red solid curve) for five cases ranging from Da≪1 to Da≫1. Here, c1=c2=1. Note that the time t is normalized by the auto-correlation time τ0 for each case. Although B(t) and A(t) share the same auto-correlation time, B(t) deviates from A(t) for cases with Da of order unity or smaller, as shown in Fig. 4a through c. On the other hand, for Da≫1, B(t) agrees with A(t) very well, as shown in Fig. 4d and e.

The simplified model has the desirable convergence property. The auto-correlation function for the simplified model (B(t) in Eq. 35) converges to that for the second version in the large-scale limit (Da→∞), as shown in Eq. (36). As confirmed in Figs. 4d and e, the two auto-correlation functions are almost identical for an integral length L greater than 10 m (or Da≥10). In the implementation of the eddy-hopping model to the LES Lagrangian cloud model, the integral length L is supposed to roughly correspond to the grid size, which is often greater than several meters to several tens of meters. Therefore, the assumption of large scales (or Da≫1) usually holds, in which case the statistical properties of the simplified model are expected to be almost unchanged after the simplification.

Figure 5 compares the time evolutions of the supersaturation fluctuations obtained from the numerical integration of the original version (dashed red curve), the second version (dotted blue curve), and the simplified model (solid green curve). Here, the numerical integration was conducted in the same manner as described in Appendix A. All results were obtained by using the same random number series. The results are shown in the time range 10τt20τ, where all cases are already in a statistically steady state.

For small scales, the simplified model produces qualitatively different trajectories of S from the second version, as shown in Fig. 5a (L=10-2 m; compare the solid green and dotted blue curves), even though these two models share the same fluctuation amplitude σS and the auto-correlation time τ0 in the statistically steady state. The difference is smaller for intermediate scales (Fig. 5b, L=100 m). For sufficiently large scales (Fig. 5c, L=102 m), the simplified model and the second version produce almost identical results.

6 Summary and conclusions

The purpose of the present paper was to obtain various statistical properties of the eddy-hopping model, a novel cloud microphysical model, which accounts for the effect of the supersaturation fluctuation at unresolved scales on the growth of cloud droplets and on spectral broadening. Two versions of the model are considered: the original version by Grabowski and Abade (2017) and the second version by Abade et al. (2018). Based on derived statistical properties, we first showed in Sect. 3 that the original version fails to reproduce a proper scaling for smaller Damköhler numbers (corresponding to small scales), resulting in a deviation of the model prediction from the reference data taken from DNSs and LESs, as shown in Fig. 1. In Sect. 4, we showed that the second version successfully reproduces the proper scaling and agrees better with the reference data than the original version for small scales (L<100 m in Fig. 2). We also showed that, by adjusting two parameters c1 and c2, the second version can almost perfectly reproduce the reference data. In Sect. 5, we discussed the possibility of simplification of the model. The simplified model consists of a single stochastic equation for the supersaturation fluctuation, as in Eq. (33), with amplitude and time parameters given by the corresponding analytical expressions for the model before the simplification. We showed that, for larger Damköhler numbers (corresponding to large scales), the auto-correlation function of the supersaturation fluctuation for the simplified model converges to that for the model before the simplification. This convergence property is desirable because the assumption of large scales usually holds in the typical parameter range for the model implementation in the LES Lagrangian cloud model.

Appendix A: Numerical integration of the eddy-hopping model

The results of the numerical integration of the original version (Eqs. 1 and 2) and the second version (Eqs. 18 and 19) are shown in Figs. 1 (blue squares) and 2 (green diamonds), respectively. For these experiments, we used the same setting as that in Sect. 5 in Thomas et al. (2020), except that the integration time was increased from 6τ to 10τ. We set a1=4.753×10-4 m−1, ε=10 cm2 s−3, and τrelax=3.513 s, and the integral time τ as

(A1) τ = 1 ( 2 π ) 1 / 3 L σ w .

As described in Appendix B, for the case of a constant dissipation rate of turbulent kinetic energy ε, σw is given as a function of L. We time integrated the governing equations of the model using 12 values of L: L=0.0128, 0.0256, 0.064, 0.128, 0.256, 0.512, 1.024, 2.56, 6.4, 12.8, 25.6, and 64.0 m. The time step δt is set as 1/1000 of τ, and the integration time is 10τ. The numerical scheme is the forward Euler method. The initial condition is such that w(0)=σwψ and S(0)=0. Each result in Figs. 1 (blue squares) and 2 (green diamonds) is obtained by averaging the results for 1000 ensembles with different seeds of random numbers.

Appendix B: Scalings for the case of a constant dissipation rate of turbulent kinetic energy

We consider classical homogeneous isotropic turbulence, in which energy is mainly injected into the system at large scales, cascaded to smaller scales by nonlinear interaction, and finally dissipated by the molecular viscosity in the smallest scales. In a statistically steady state, the dissipation rate of turbulent kinetic energy is defined as ε. If ε is fixed and the integral scale L is changed, then the kinetic energy E scales as follows (Thomas et al.2020):

(B1) E ( L ε ) 2 / 3 .

The black dots in Fig. B1 show the relation between L and E in the reference data taken from DNSs and LESs by Thomas et al. (2020) (Table 2 of their study). In their simulation, the dissipation rate was fixed to ε=10 cm2s−3. The orange curve in Fig. B1 indicates the function E=αε2/3L2/3, where α is the fitting parameter. The best fit is given by α=0.475. The root-mean-square turbulent velocity is calculated as a function of L by urms=σw=(2E/3), and σw is used as the parameter in the eddy-hopping model. Note that Thomas et al. (2020) used the same type of large-scale forcing as that used by Kumar et al. (2012), where the integral length L is set to be equal to the box length Lbox.

https://acp.copernicus.org/articles/21/13119/2021/acp-21-13119-2021-f06

Figure B1 Relationship between the integral scale L and the turbulent kinetic energy E. The black dots are taken from the reference data in Thomas et al. (2020). The orange curve indicates the fitting function E=αε2/3L2/3 with α=0.475.

Download

Appendix C: Achievement of a statistically steady state

We confirm that all of the results of the numerical integration of the eddy-hopping model in the present study achieved statistically steady states. For this purpose, we first derive the analytical expression for the time evolutions of the variance and covariance of the variables in the model and then compare these analytical expressions with the results of the numerical integration.

The governing equations given by Eqs. (3) and (2) can be rewritten in generalized forms as

(C1)dwdt=-1τ1w(t)+Fw(t),(C2)dSdt=a1w(t)-S(t)τ2,

where τ1 and τ2 are the relaxation times for w and S, respectively, and the forcing term Fw(t) satisfies Eq. (4). Evolution equations for the variance and covariance of the variables are derived as follows:

(C3)dVw(t)dt=-2τ1Vw(t)+2σw2τ1,(C4)dC(t)dt=a1Vw(t)-1τ1+1τ2C(t),(C5)dVS(t)dt=-2τ2VS(t)+2a1C(t),

where Vw(t), C(t), and VS(t) are, respectively, defined as

(C6)Vw(t)=w(t)w(t),(C7)C(t)=w(t)S(t),(C8)VS(t)=S(t)S(t).

For the numerical integration of the eddy-hopping model by Thomas et al. (2020), τ1=τ and τ2=τrelax. Since the initial conditions for w(t) and S(t) are set to w(0)=σwψ and S(0)=0 in Thomas et al. (2020), the corresponding initial conditions for the variance and covariance are given by

(C9)Vw(0)=σw2,(C10)C(0)=0,(C11)VS(0)=0.

Solving Eqs. (C3) through (C5) with the initial conditions given by Eqs. (C9) through (C11), we obtain

(C12)Vw(t)=σw2,(C13)C(t)=a1σw2τ31-e-t/τ3,(C14)VS(t)=a12σw2τ3τ21-e-2t/τ2+2a12σw2τ3τ4e-t/τ3-e-2t/τ2,

where τ3 and τ4 are, respectively, defined as

(C15) τ 3 = τ 1 τ 2 τ 1 + τ 2 , and τ 4 = τ 1 τ 2 τ 2 - τ 1 .
https://acp.copernicus.org/articles/21/13119/2021/acp-21-13119-2021-f07

Figure C1 Standard deviation of the supersaturation fluctuation σS at times t=0.6τ (a) and t=6τ (b) obtained from the analytical expression given by Eq. (C14) (cyan dots) and the results of the numerical integration of the original version given by Eqs. (1) and (2) (black crosses). The orange curve, red triangles, and axes of the panel are the same as in Fig. 1. The two short black lines indicate slopes of 2/3 and 1/3. The setting for the numerical integration is the same as that used in Sect. 3, except that the integration times are 0.6τ and 6τ in (a) and (b), respectively.

Download

Figure C1 compares the analytical expression given by Eq. (C14) (cyan dots) with the results of the numerical integration of the original version given by Eqs. (1) and (2) (black crosses). The setting for the numerical experiment is the same as that used in Fig. 1, except that the integration time is 0.6τ in Fig. C1a and 6τ in Fig. C1b. The results of the numerical integration (black crosses) agree well with the analytical expression (cyan dots), and both approach the theoretical curve for the statistically steady state (orange curve in each panel) as the integration time increases. Figure C1a also indicates that the results of the numerical integration of the eddy-hopping model by Thomas et al. (2020) (red triangles) are fairly close to our results at 0.6τ. Thus, it might be possible that the integration time of their numerical experiment was not long enough to achieve a statistically steady state.

Appendix D: Derivation of auto-correlation function

We derive the analytical expression for the auto-correlation function of the supersaturation fluctuation S(t) in the eddy-hopping model. As in Appendix C, we start from the generalized form of the eddy-hopping model as follows:

(D1)dwdt=-1τ1w(t)+Fw(t),(D2)dSdt=a1w(t)-S(t)τ2,

where τ1 and τ2 are the relaxation times for w and S, respectively, and the forcing term Fw(t) satisfies Eq. (4). We consider that the system is in a statistically steady state.

First, multiplying Eq. (D2) by et/τ2 and applying the product rule of differentiation, we obtain

(D3) d d t S ( t ) e t / τ 2 = a 1 w ( t ) e t / τ 2 .

Integrating Eq. (D3) from t=0 to t, we obtain

(D4) S ( t ) = S ( 0 ) e - t / τ 2 + 0 t a 1 w ( ξ ) e ( ξ - t ) / τ 2 d ξ .

(Note that we chose the integration range [0,t] for simplicity of notation. Since we consider a statistically steady state, the following discussion is unchanged if the integration range is [t0,t0+t].) Applying a similar procedure to that above to Eq. (D1) with the integration range t:0ξ, we obtain

(D5) w ( ξ ) = w ( 0 ) e - ξ / τ 1 + 0 ξ F w ( ζ ) e ( ζ - ξ ) / τ 1 d ζ .

Substituting Eq. (D5) into (D4) and calculating some of the integrations, we obtain

(D6)S(t)=S(0)e-t/τ2+0ta1(w(0)e-ξ/τ1+0ξFw(ζ)e(ζ-ξ)/τ1dζ)e(ξ-t)/τ2dξ(D7)=S(0)e-t/τ2+a1w(0)e-t/τ20te(τ2-1-τ1-1)ξdξ+a1e-t/τ20t0ξFw(ζ)eζ/τ1e(τ2-1-τ1-1)ξdζdξ(D8)=S(0)e-t/τ2+a1w(0)τ2-1-τ1-1-1e-t/τ1-e-t/τ2+a1e-t/τ20t0ξFw(ζ)eζ/τ1e(τ2-1-τ1-1)ξdζdξ.

Multiplying Eq. (D8) by S(0) and taking an ensemble average, we obtain

(D9) S ( t ) S ( 0 ) = S ( 0 ) S ( 0 ) e - t / τ 2 + a 1 w ( 0 ) S ( 0 ) τ 2 - 1 - τ 1 - 1 - 1 e - t / τ 1 - e - t / τ 2 ,

because of the statistical independence Fw(ζ)S(0)=0. Next, as in the derivation of Eq. (11), we multiply Eq. (D2) by S and consider the statistically steady state. We obtain

(D10) S ( 0 ) S ( 0 ) = a 1 τ 2 w ( 0 ) S ( 0 ) .

Substituting Eq. (D10) into (D9), we have

(D11) S ( t ) S ( 0 ) = S ( 0 ) S ( 0 ) e - t / τ 2 + S ( 0 ) S ( 0 ) τ 2 - 1 τ 2 - 1 - τ 1 - 1 - 1 e - t / τ 1 - e - t / τ 2 .

Therefore, the auto-correlation function of the supersaturation fluctuation S(t) for the eddy-hopping model given by Eqs. (D1) and (D2) in the statistically steady state is written as follows:

(D12)A(t)=S(t)S(0)S(0)S(0)(D13)=e-t/τ2+τ1τ1-τ2e-t/τ1-e-t/τ2(D14)=τ1τ1-τ2e-t/τ1-τ2τ1-τ2e-t/τ2.

The auto-correlation time τ0 is obtained by time-integrating A(t) as

(D15)τ0=0A(t)dt=τ12-τ22τ1-τ2(D16)=τ1+τ2.

For the original version of the eddy-hopping model given by Eqs. (1) and (2), we have

(D17) τ 1 = τ , and τ 2 = τ relax .

For the second version given by Eqs. (18) and (19), we have

(D18) τ 1 = c 1 τ , and τ 2 = 1 c 1 τ + 1 c 2 τ relax - 1 .
Code availability

Programs for numerical integration of the eddy-hopping model are written in Fortran 90 and are available upon request.

Data availability

The data in this study are available upon request.

Author contributions

IS conducted the numerical simulations and data analysis. IS and TW performed the theoretical analyses in Sect. 3 and 4. IS and TG performed the theoretical analysis in Sect. 5. All three of the authors were involved in preparing the manuscript.

Competing interests

The authors declare that they have no conflicts of interest.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

We are grateful to Kei Nakajima for his technical support. The present study was supported by MEXT KAKENHI grant nos. 20H00225 and 20H02066, by JSPS KAKENHI grant no. 18K03925, by the Naito Foundation, by the HPCI System Research Project (project ID: hp200072, hp210056), by the NIFS Collaboration Research Program (NIFS20KNSS143), by the Japan High Performance Computing and Networking plus Large-scale Data Analyzing and Information Systems (jh200006, jh210014), and by High Performance Computing (HPC 2020, HPC2021) at Nagoya University.

Review statement

This paper was edited by Markus Petters and reviewed by three anonymous referees.

References

Abade, G. C., Grabowski, W. W., and Pawlowska, H.: Broadening of cloud droplet spectra through eddy hopping: Turbulent entraining parcel simulations, J. Atmos. Sci., 75, 3365–3379, https://doi.org/10.1175/JAS-D-18-0078.1, 2018. a, b, c, d, e, f, g, h, i

Chandrakar, K. K., Cantrell, W., Chang, K., Ciochetto, D., Niedermeier, D., Ovchinnikov, M., Shaw, R. A., and Yang, F.: Aerosol indirect effect from turbulence-induced broadening of cloud-droplet size distributions, P. Natl. Acad. Sci. USA, 113, 14243–14248, https://doi.org/10.1073/pnas.1612686113, 2016. a

Chandrakar, K. K., Grabowski, W. W., Morrison, H., and Bryan, G. H.: Impact of entrainment–mixing and turbulent fluctuations on droplet size distributions in a cumulus cloud: An investigation using Lagrangian microphysics with a sub–grid–scale model, J. Atmos. Sci., https://doi.org/10.1175/JAS-D-20-0281.1, 2021. a

Clark, T. L. and Hall, W. D.: A numerical experiment on stochastic condensation theory, J. Atmos. Sci., 36, 470–483, https://doi.org/10.1175/1520-0469(1979)036<0470:ANEOSC>2.0.CO;2, 1979. a

Cooper, W. A.: Effects of variable droplet growth histories on Droplet size distributions. Part I: Theory, J. Atmos. Sci., 46, 1301–1311, https://doi.org/10.1175/1520-0469(1989)046<1301:EOVDGH>2.0.CO;2, 1989. a

Devenish, B. J., Bartello, P., Brenguier, J.-L., Collins, L. R., Grabowski, W. W., IJzermans, R. H. A., Malinowski, S. P., Reeks, M. W., Vassilicos, J. C., Wang, L.-P., and Warhaft, Z.: Droplet growth in warm turbulent clouds, Q. J. Roy. Meteor. Soc., 138, 1401–1429, https://doi.org/10.1002/qj.1897, 2012. a

Grabowski, W. W. and Abade, G. C.: Broadening of cloud droplet spectra through eddy hopping: Turbulent adiabatic parcel simulations, J. Atmos. Sci., 74, 1485–1493, https://doi.org/10.1175/JAS-D-17-0043.1, 2017. a, b, c, d, e, f, g

Grabowski, W. W. and Wang, L.-P.: Growth of cloud droplets in a turbulent environment, Annu. Rev. Fluid Mech., 45, 293–324, https://doi.org/10.1146/annurev-fluid-011212-140750, 2013. a

Korolev, A. V. and Mazin, I. P.: Supersaturation of water vapor in clouds, J. Atmos. Sci., 60, 2957–2974, 2003. a, b

Kostinski, A. B.: Simple approximations for condensational growth, Environ. Res. Lett., 4, 015005, https://doi.org/10.1088/1748-9326/4/1/015005, 2009. a

Kumar, B., Janetzko, F., Schumacher, J., and Shaw, R. A.: Extreme responses of a coupled scalar–particle system during turbulent mixing, New J. Phys., 14, 115020, https://doi.org/10.1088/1367-2630/14/11/115020, 2012. a

Lanotte, A. S., Seminara, A., and Toschi, F.: Cloud droplet growth by condensation in homogeneous isotropic turbulence, J. Atmos. Sci., 66, 1685–1697, https://doi.org/10.1175/2008JAS2864.1, 2009. a, b

Lasher-Trapp, S. G., Copper, W. A., and Blyth, A. M.: Broadening of droplet size distributions from entrainment and mixing in a cumulus cloud, Q. J. Roy. Meteor. Soc., 131, 195–220, https://doi.org/10.1256/qj.03.199, 2005. a

Marcq, P. and Naert, A.: A Langevin equation for the energy cascade in fully developed turbulence, Physica D, 124, 368–381, https://doi.org/10.1016/S0167-2789(98)00237-1, 1998. a

Paoli, R. and Shariff, K.: Turbulent condensation of droplets: Direct simulation and a stochastic model, J. Atmos. Sci., 66, 723–740, https://doi.org/10.1175/2008JAS2734.1, 2009. a

Politovich, M. K. and Cooper, W. A.: Variability of the supersaturation in cumulus clouds, J. Atmos. Sci., 45, 2064–2086, https://doi.org/10.1175/1520-0469(1988)045<1651:VOTSIC>2.0.CO;2, 1988. a

Pope, S. B.: Turbulent flows, Cambridge University Press, Cambridge, 2000. a

Saito, I., Gotoh, T., and Watanabe, T.: Broadening of cloud droplet size distributions by condensation in turbulence, J. Meteor. Soc. Japan, 97, 867–891, https://doi.org/10.2151/jmsj.2019-049, 2019a. a

Saito, I., Watanabe, T., and Gotoh, T.: A new time scale for turbulence modulation by particles, J. Fluid Mech., 880, R6, https://doi.org/10.1017/jfm.2019.775, 2019b. a

Sardina, G., Picano, F., Brandt, L., and Caballero, R.: Continuous growth of droplet size variance due to condensation in turbulent clouds, Phys. Rev. Lett., 115, 1–5, https://doi.org/10.1103/PhysRevLett.115.184501, 2015. a, b

Sardina, G., Picano, F., Brandt, L., and Caballero, R.: Broadening of cloud droplet size spectra by stochastic condensation: effects of mean updraft velocity and CCN activation, J. Atmos. Sci., 75, 451–467, https://doi.org/10.1175/JAS-D-17-0241.1, 2018. a

Sawford, B. L.: Reynolds number effects in Lagrangian stochastic models of turbulent dispersion, Phys. Fluids A-Fluid, 3, 1577–1586, https://doi.org/10.1063/1.857937, 1991. a

Sedunov, Y. S.: Physics of Drop Formation in the Atmosphere, John Wiley and Sons, New York, 1974.  a

Shaw, R. A.: Particle-turbulence interactions in atmospheric clouds, Annu. Rev. Fluid Mech., 35, 183–227, https://doi.org/10.1146/annurev.fluid.35.101101.161125, 2003. a

Siewert, R., Bec, J., and Krstulvic, G.: Statistical steady state in turbulent droplet condensation, J. Fluid Mech., 810, 254–280, https://doi.org/10.1017/jfm.2016.712, 2017. a

Thomas, L., Grabowski, W. W., and Kumar, B.: Diffusional growth of cloud droplets in homogeneous isotropic turbulence: DNS, scaled-up DNS, and stochastic model, Atmos. Chem. Phys., 20, 9087–9100, https://doi.org/10.5194/acp-20-9087-2020, 2020. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p

Download
Short summary
We provide various statistical properties for the stochastic model of eddy hopping, which is a novel cloud microphysical model that accounts for the effect of the supersaturation fluctuation at unresolved scales on the growth of cloud droplets and on spectral broadening in a turbulent cloud. Our results indicate that the model can be improved to have better fidelity to the reference data and to require less computational cost.
Altmetrics
Final-revised paper
Preprint