Statistical properties of a stochastic model of eddy hopping

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. It is shown that the model fails to reproduce a proper scaling for a certain range of parameters, resulting in a deviation of the model prediction from the reference data taken from direct numerical simulations and large-eddy simulations 5 (LESs). Corrections to the model are introduced so that the corrected model can accurately reproduce the reference data with the proper scaling. In addition, a possible simplification of the model is discussed, which may contribute to a reduction in computational cost while keeping the statistical properties almost unchanged in the typical parameter range for the model implementation in the LES Lagrangian cloud model.

the possibility of simplification of the model, which may contribute to a reduction in computational cost 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. Section 3 presents 30 a theoretical analysis and numerical experiments and demonstrates the improper scaling in the model prediction. Section 4 introduces corrections to the model. Finally, Section 5 discusses the possibility of simplification of the model.

Governing equations
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 35 process: 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 Here, the first term on the right-hand side represents the effect of adiabatic cooling/warming due to air parcel ascent/descent caused by the vertical velocity w ′ (t). The parameter a 1 has the unit of a scalar gradient. The second term on the right-hand side represents the effect of condensation/evaporation of droplets. The time scale τ 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 Cooper, 1988; 45 Korolev and Mazin, 2003;Kostinski, 2009;Devenish et al., 2012).
Equation (1) can also be written as the following derivative form (Pope, 2000): Here, the term F w ′ (t) is statistically independent of S ′ and obeys the Gaussian random process that has zero mean and two-time 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 eddy-hopping model.
We now obtain the analytical expression for the standard deviation of the supersaturation fluctuation, σ S ′ , in a statistically 55 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 because of statistical independence (⟨S ′ F w ′ ⟩ = 0). Second, multiplying Eq. (2) by w ′ and taking an ensemble average, we Summing Eqs. (5) and (6), we obtain Next, we consider a statistically steady state. Since an ensemble-averaged variable does not change in time (d ⟨•⟩ /dt = 0) and ⟨w ′2 ⟩ = σ 2 w ′ in the statistically steady state, we obtain the flux of the supersaturation in the vertical direction as follows: where Da is the Damköhler number (Shaw, 2003) defined as Next, multiplying Eq. (2) by S ′ and taking an ensemble average, we obtain In the statistically steady state, we have Combining Eqs. (8) and (11), we obtain or equivalently, Here, σ S ′ in Eq. (13) has two important asymptotic forms, as shown below: For the case of a constant dissipation rate of turbulent kinetic energy ε, σ w ′ ∼ L 1/3 (see Appendix B), and we have the 85 following scaling:

90
For the case of a constant dissipation rate of turbulent kinetic energy ε, σ w ′ ∼ L 1/3 and τ ∼ L 2/3 (see Appendix B), and we have the following scaling: 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 σ S ′ ∼ a 1 τ relax σ w ′ for the large scale limit, which is consistent 95 with Eq. (14). On the other hand, we should have σ S ′ ∼ a 1 τ σ w ′ for the small scale limit, which is inconsistent with Eq. (16).
Therefore, the eddy-hopping model given by Eqs. (3) and (2) does not reproduce the proper scaling for the small scale limit. 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 eddy-hopping model given by Eqs.
(1) and (2) (blue squares). Here, numerical integration is conducted in the same manner as that by Thomas et al. (2020) (Section 5 of their study), except that the integration time 100 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 σ S ′ ∼ L 1/3 for large scales (approximately L > 10 1 m) and the improper scaling σ S ′ ∼ L 2/3 for small scales (approximately L < 10 −1 m). These results are contrary to the results of DNSs and LESs (scaled-up DNSs) conducted by Thomas et al. (2020) (black circles in Figure 1), which show the proper 105 scalings both for large and small scales (σ S ′ ∼ L 1/3 and ∼ L 1 , respectively).
Note that Figure 1 also shows the results of the numerical integration of the eddy-hopping model 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. is the integral length L. The black circles 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 eddy-hopping model reported by Thomas et al. (2020). The range of L and σ S ′ for the panel is the same as in Figure 10 in Thomas et al. (2020). The three short black lines indicate slopes of 1, 2/3, and 1/3.
The eddy-hopping model given by Eqs. (1) and (2) shows the improper scaling for small scales because of the assumption 110 made in the formulation of the model. Originally, Eq.
(2) [corresponding to Eq. (8) in Grabowski and Abade (2017)] 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 eddy-hopping model does not reproduce the proper scaling for small scales.
In the following section, we introduce some corrections to the original model given by Eqs. (1) and (2) so that the corrected model can reproduce proper scalings both for large and small scales. 115 We introduce corrections to the eddy-hopping model so that the model can reproduce the reference data taken from DNSs and LESs in Figure 1. We first show the results. Equations (1) and (2) are, respectively, corrected as follows: There are two types of corrections. First, we added a term that is proportional to −S ′ /τ to 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(Sardina et al., , 2018Chandrakar et al., 2016;Siewert et al., 2017;Saito et al., 2019a). The time scale of the damping effect due to turbulent mixing is characterized by the integral time τ , whereas that due to condensation/evaporation of cloud droplets is characterized by the phase relaxation time τ relax . The 125 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). We also note that the term −S ′ /τ has been introduced to the equation for S ′ in the eddy-hopping model by Abade et al. (2018) [see Eq.
(15) of their study]. However, a subsequent study by Thomas et al. (2020) used a model without the term −S ′ /τ [see Eq. (14) of their study]. Although there might be some confusion, we again emphasize the importance of this term, which plays an essential role in reproducing the 130 reference data, as described below.
The second correction to the model is that we changed two time scales, τ and τ relax , to (c 1 τ ) and (c 2 τ relax ), respectively.
Here, c 1 and c 2 are constants and are used as tuning parameters. These types of parameters are not new. For example, a parameter corresponding to c 1 is commonly used in the Langevin stochastic equation in turbulence research (Sawford, 1991;Marcq and Naert, 1998). Formally, the inverse of c 1 is referred to as the drift coefficient, and the coefficients for the velocity 135 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 c 2 has been demonstrated in a recent study on turbulence modulation by particles (Saito et al., 2019b) Applying the analytical procedure described in Section 3 to the corrected model given by Eqs. (18) and (19), we first obtain instead of Eq. (8). Next, instead of Eq. (11), we have ( 1 Finally, the analytical expression corresponding to Eq. (13) is Asymptotic forms of σ S ′ in Eq. (22) for the large and small scale limits are, respectively, given as follows: For the case of a constant dissipation rate of turbulent kinetic energy ε, we have 2. Small scale limit which indicates that τ 1/2 relax 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 Therefore, the corrected model reproduces asymptotic forms σ S ′ ∼ a 1 τ relax σ w ′ and ∼ a 1 τ σ w ′ for the large and small scale 160 limits, respectively, which are both consistent with the result of the scaling argument by Lanotte et al. (2009).
Two parameters c 1 and c 2 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 c 1 = 0.746 and c 2 = 1.28. Figure 2 (orange curve) shows the theoretical curve given by Eq. (22) with these values of c 1 and c 2 , which agrees well with the reference data (black circles). Figure 2 (green diamonds) also shows the results of the numerical integration of the corrected model, which agree with the 165 theoretical curve, as expected. Here, the numerical integration was conducted in the same manner as in the previous section (see Appendix A for details). These results confirm that the corrected eddy-hopping model given by Eqs. (18) and (19), with proper values of c 1 and c 2 , can accurately reproduce the reference data.

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 corrected 170 model given by Eqs. (18) and (19), but the same argument also applies to the original model given by Eqs. (1) and (2).
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 Section 4, we obtained an analytical 175 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:

180
where τ 1 and τ 2 are, respectively, defined as 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 The auto-correlation function A(t) in Eq. (28) respectively. Second, for the small scale limit (Da → 0), the asymptotic forms of A(t) and τ 0 are given by e −t/(c1τ ) and lim respectively.
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:

195
where σ S ′ and τ 0 are given by Eqs. (22) and (30), respectively. Note that the simplified model given by Eq. (33) is a singleequation model, as compared to the two-equation model given by Eqs. (18) and (19) before the simplification. The autocorrelation function for S ′ in the simplified model given by Eq. (33) is given by 200 which has the following two asymptotic forms. First, for the large scale limit (Da → ∞), which agrees with the corresponding asymptotic form given by Eq. (31) for the corrected model. Second, for the small scale   3(e), the two auto-correlation functions are almost identical for an integral length L greater than 10 m (or Da ≥ 10). In 220 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.

Summary and conclusions 225
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. Based on derived statistical properties, we first showed in Section 3 that the model 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 Figure 1. In Section 4, we introduced two 230 corrections to the model so that the corrected can accurately reproduce the reference data with the proper scaling. The first correction is to add the term representing the effect of turbulent mixing (eddy diffusivity) on the supersaturation fluctuation, and the second is to multiply time scales τ (integral time) and τ relax (phase relaxation time) by tuning parameters c 1 and c 2 , respectively. In Section 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 235 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. Since the assumption of large scales usually holds in the typical parameter range for the model implementation in the LES Lagrangian cloud model, the simplified model may contribute to a reduction in computational cost while keeping the statistical properties almost unchanged after the simplification.

240
In a future study, the actual performance of the simplified model should be validated through numerical simulations of cloud models, such as those by Grabowski and Abade (2017); Abade et al. (2018).

Appendix A: Numerical integration of the eddy-hopping model
The results of the numerical integration of the eddy-hopping model [Eqs.
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, 250 0.512, 1.024, 2.56, 6.4, 12.8, 25.6, and 64.0 m. The time step δt is set as 1/1,000 of τ , and the integration time is 10τ . In addition, the numerical scheme is the forward Euler method. The initial condition is such that w ′ (0) = σ w ′ ψ and S ′ (0) = 0.
Each result in Figures 1 (blue squares) and 2 (green diamonds) is obtained by averaging the results for 1,000 ensembles with different seeds of random numbers.
Appendix B: Scalings for the case of a constant dissipation rate of turbulent kinetic energy 255 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): The black circles in Figure 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 cm 2 s −3 . The orange curve in Figure B1 indicates the function E = αε 2/3 L 2/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 u rms = σ 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 265 Kumar et al. (2012), where the integral length L is set to be equal to the box length L box . Fit DNS+LES TG20 Figure B1. Relationship between the integral scale L and the turbulent kinetic energy E. The black circles are taken from the reference data in Thomas et al. (2020). The orange curve indicates the fitting function E = αε 2/3 L 2/3 with α = 0.475.

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 inte-270 gration.
The governing equations given by Eqs. (3) and (2) can be rewritten in generalized forms as where τ 1 and τ 2 are the relaxation times for w ′ and S ′ , respectively, and the forcing term F w ′ (t) satisfies Eq. (4). Evolution 275 equations for the variance and covariance of the variables are derived as follows: where V w ′ (t), C(t), and V S ′ (t) are, respectively, defined as 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 285 for the variance and covariance are given by Solving Eqs. (C3) through (C5) with the initial conditions given by Eqs. (C9) through (C11), we obtain where τ 3 and τ 4 are, respectively, defined as Figure C1 compares the analytical expression given by Eq. (C14) (cyan circles) with the results of the numerical integration of the eddy-hopping model given by Eqs.
(1) and (2) (black crosses). The setting for the numerical experiment is the same as that used in Figure 1, except that the integration time is 0.6τ in Figure C1(a) and 6τ in Figure C1(b). The results of the numerical integration (black crosses) agree well with the analytical expression (cyan circles), and both approach the theoretical curve for the statistically steady state (orange curve in each panel) as the integration time increases. Figure C1(a) also indicates 300 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. (1) and (2) (black crosses) The orange curve, red triangles, and axes of the panel are the same as in Figure 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 Section 3, except that the integration times are 0.6τ and 6τ in (a) and (b), respectively.

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-305 hopping model. As in Appendix C, we start from the generalized form of the eddy-hopping model as follows: where τ 1 and τ 2 are the relaxation times for w ′ and S ′ , respectively, and the forcing term F w ′ (t) satisfies Eq. (4). We consider that the system is in a statistically steady state.
Substituting Eq. (D5) into Eq. (D4) and calculating some of the integrations, we obtain Multiplying Eq. (D8) by S ′ (0) and taking an ensemble average, we obtain because of the statistical independence (⟨F w ′ (ζ)S ′ (0)⟩ = 0). Next, as in the derivation of Eq. (11) and (D2) in the statistically steady state is written as follows: The auto-correlation time τ 0 is obtained by time-integrating A(t) as For the original version of the eddy-hopping model given by Eqs.