ABSTRACT: Turbulent fluctuations of scalar and velocity fields are critical for cloud microphysical processes, e.g., droplet activation and size distribution evolution, and can therefore influence cloud radiative forcing and precipitation formation. Lagrangian and Eulerian water vapor, temperature, and supersaturation statistics are investigated in direct-numerical simulations (DNS) of turbulent Rayleigh-Bénard convection in the Pi Convection Cloud Chamber to provide a foundation for parameterizing subgrid-scale fluctuations in atmospheric models. A subgrid model for water vapor and temperature variances and covariance and supersaturation variance is proposed, valid for both clear and cloudy conditions. Evaluation of phase change contributions through an a priori test using DNS data shows good performance of the model. Supersaturation is a nonlinear function of temperature and water vapor, and relative external fluxes of water vapor and heat (e.g., during entrainment-mixing and phase change) influence turbulent supersaturation fluctuations. Although supersaturation has autocorrelation and structure functions similar to the independent scalars (temperature and water vapor), the autocorrelation timescale of supersaturation differs. Relative scalar fluxes in DNS without cloud make supersaturation PDFs less skewed than the adiabatic case, where they are highly negatively skewed. However, droplet condensation changes the PDF shape response: it becomes positively skewed for the adiabatic case and negatively skewed when the sidewall relative fluxes are large. Condensation also increases correlations between water vapor and temperature in the presence of relative scalar fluxes but decreases correlations for the adiabatic case. These changes in correlation suppress supersaturation variability for the non-adiabatic cases and increase it for the adiabatic case. Implications of this work for subgrid microphysics modeling using a Lagrangian stochastic scheme are also discussed.

Turbulent fluctuations in scalar and dynamical fields are critical for cloud microphysical processes. For example, supersaturation and the associated variability due to turbulent stirring and mixing is crucial for condensation growth of droplets. Past studies have highlighted the importance of supersaturation variability in droplet size distribution broadening and droplet activation (e.g., Mazin 1968;Cooper 1989;Korolev 1995;Paoli and Shariff 2009;Sardina et al. 2015;Chandrakar et al. 2016Chandrakar et al. , 2017;;Siewert et al. 2017;Grabowski and Abade 2017;Saito et al. 2019;Prabhakaran et al. 2020). These impacts of supersaturation variability on droplet size distribution can potentially influence the macroscopic cloud properties and processes: cloud radiative forcing, precipitation formation, and feedback to the dynamics. Supersaturation depends linearly on water vapor and non-linearly on temperature. Thus, treating it like a single independent scalar (e.g., in Lanotte et al. 2009;Sardina et al. 2015;Siewert et al. 2017;Gotoh et al. 2021) could lead to errors in its evolution and interaction with cloud droplets.

Even in the absence of external fluxes, the different diffusivities of water vapor and temperature affect the statistics of supersaturation variability (Chandrakar et al. 2020;Thomas et al. 2021).

Moreover, the presence of external water vapor and temperature fluxes, for example, from turbulent entrainment of environmental air into clouds, affects supersaturation fluctuations differently from those of independent scalars (Wyngaard et al. 1978;Siebert and Shaw 2017;Chandrakar et al. 2022). Chandrakar et al. (2022) showed that these relative scalar fluxes decorrelate water vapor and temperature fields and increase supersaturation fluctuations in moist convection without cloud.

The current study further extends this investigation to cloudy conditions. How can we account for these turbulence effects on clouds in models across scales? One option is a bottom-up approach where information from controlled experiments/observations or finescale direct-numerical simulations (DNS) provides a foundation for subgrid modeling in coarserresolution models (e.g., cloud-resolving or large-eddy-simulation (LES)). Chandrakar et al. (2022) introduced a new subgrid supersaturation model based on scalar variance and covariance prognosis (see Appendix A here) and evaluated it against the benchmark DNS of the Pi Convection Cloud Chamber (Chang et al. 2016) with sidewall fluxes mimicking entrainment. The proposed subgrid model performs exceptionally well in an LES framework to reproduce DNS results. However, the model only considered the case without clouds. The current study extends this model to cloudy conditions. It accounts for the phase change contributions (discussed in Sec. 2) to water vapor and temperature variances and their covariance. Moreover, this study investigates assumptions in subgrid models of turbulence-microphysics interactions for LES and cloud-resolving models using DNS of turbulent Rayleigh-Bénard convection with droplets in the Pi Convection Cloud Chamber.

Rayleigh-Bénard convection in the chamber, which is well-studied and simpler than atmospheric clouds with complex feedbacks, is well suited for model development. Moreover, well-controlled turbulence and cloud droplet measurements in the chamber provide a reference for simulation studies.

The Lagrangian stochastic approach based on the Langevin drift-diffusion equation is one of the methods employed in the past to couple Eulerian subgrid models (e.g., the one discussed above) with microphysical schemes and represent interactions of cloud particles with fluctuating scalar fields (water vapor and temperature) in turbulent clouds. It assumes scalar fluctuations as a Markov process where the future state depends only on the present state, not the past. It broadly has two components: a deterministic drift component that relaxes the fluctuations toward the reference state with a characteristic timescale of a turbulent field (𝜏 𝑡 ); and the diffusion driver (Wiener forcing) that represents the unresolved turbulent gradient of a scalar as a delta-correlated Gaussian random process. A simple example of this model for a scalar (𝜙) in a turbulent flow is shown below:

𝜏 𝑡 𝑑𝑡 + 𝐴(𝑥, 𝑡)𝑑𝑊 (𝑡).

(1)

where 𝑥 and 𝑡 are the location and time and prime denotes fluctuation. The drift term (the first term in the above equation) ensures an exponential autocorrelation function while turbulent mixing (a criterion for the Markov process) and the increment of the Wiener process (the second term)

drives small-scale fluctuations with the assumption of a linear second order-structure function and magnitude 𝐴(𝑥, 𝑡) at location 𝑥 and time 𝑡. The Wiener process 𝑊 is a continuous random Gaussian process with time-independent Gaussian increments having zero mean and finite variance. Please refer to the Gardiner (2004) textbook for a comprehensive overview of the Langevin model and stochastic differential equations. The magnitude of the Wiener forcing term (𝐴(𝑥, 𝑡)) at the particle location can be obtained from the Eulerian subgrid model introduced in Chandrakar et al. (2022) and given here in Appendix A as well.

In this study, assumptions involved in Lagrangian stochastic models of supersaturation variability are evaluated. This work provides a foundation for future improvement of such models. Common assumptions in studies of stochastic condensation and its subgrid representation in models are:

• Supersaturation variability behaves similarly to independent scalars (water vapor and temperature) (e.g., Lanotte et al. 2009;Sardina et al. 2015;Siewert et al. 2017;Gotoh et al. 2021).

• The PDF of supersaturation fluctuations (or Lagrangian temporal differences) follows a Gaussian shape.

• In the Lagrangian stochastic model, forcing (Wiener forcing) of supersaturation fluctuations, or Lagrangian temporal differences, is delta correlated in time.

• The Lagrangian autocorrelation function of supersaturation is exponential (a criterion for the Markov process), similar to the velocity and independent scalar fields. (Tennekes 1979) • The Lagrangian integral timescale of supersaturation is the same as the integral timescale of water vapor, temperature, and velocity fields.

• The second-order Lagrangian structure function of supersaturation follows the Lagrangian equivalent of the Kolmogorov-Obukhov-Corrsin scaling: ⟨𝑑𝜙 2 (𝜏)⟩ = 𝐶𝜀 𝜙 𝜏 (Monin and Yaglom 1971;Tennekes 1979;Biferale et al. 2008). Here 𝜙 is a scalar/velocity field, 𝜏 is the time lag, 𝐶 is a universal constant, and 𝜀 is the dissipation rate.

The shape of the supersaturation PDF is critical for subgrid modeling of drop condensation. Gaussian mixing model. The supersaturation PDF becomes more symmetric when the diffusivity difference between water vapor and temperature is artificially increased (Chandrakar et al. 2020).

Further investigation of supersaturation PDFs using DNS is required where smaller scales can be appropriately resolved, especially in the presence of relative external fluxes of water vapor and heat. Moreover, the presence of cloud particles introduces relative fluxes of scalars through phase changes in turbulent moist convection. Latent heating adds thermal energy to the system while condensation removes water vapor. This introduces additional sources/sinks to water vapor and temperature variability and also affects their correlation. Therefore, it is critical to understand the influence of phase changes on scalar and supersaturation statistics for appropriately modeling the cloud-turbulence interactions. This article further expands on this topic and presents Eulerian and Lagrangian supersaturation statistics from DNS and their comparison with independent scalar statistics, and also investigates the impact of phase change via condensation and evaporation.

Another aspect of droplet-turbulence interaction is droplet inertia and the presence of gravity that lead to deviation in particle trajectories in turbulent flow compared to Lagrangian passive parcels. The impact of inertia for small cloud droplets (not rain and drizzle drops) is typically ignored for most cloud conditions since the turbulent Stokes number (ratio of droplet inertial relaxation timescale to Kolmogorov timescale of small eddies) is generally ≪ 0.1. However, the impact of gravitational settling of cloud droplets could be significant, especially at small scales.

The presence of gravitational settling could act as a filter for small-scale eddies, decorrelating the supersaturation field faster along droplet trajectories and affecting the covariance of water vapor and temperature. This is a point of concern for models of subgrid cloud-turbulence interactions, and it is typically ignored. In this study, we address the topic by analyzing Lagrangian trajectories of passive (without mass and thermal inertia) and inertial (with mass but without thermal inertia) cloud particles in the DNS.

Following the above discussion, in this article we focus on the following science questions:

• Does supersaturation display characteristics similar to independent scalars in turbulent convection?

• Characteristics of Lagrangian turbulence statistics (e.g., Lagrangian autocorrelation function, second-order structure-function, and PDFs of supersaturation gradient/differences) are critical for subgrid-scale modeling. How do these Lagrangian statistics of supersaturation fluctuations behave relative to water vapor and temperature fluctuations in turbulent cloudy convection?

• How does cloud droplet settling affect supersaturation correlations and Lagrangian statistics?

• How does the shape of the supersaturation fluctuation PDF change across different relative flux conditions? Is the Gaussian shape of the supersaturation PDF assumed in previous subgrid models (e.g., Paoli and Shariff 2009;Sardina et al. 2015;Chandrakar et al. 2016;Grabowski and Abade 2017) a good approximation under realistic conditions?

• How does the condensation process affect scalar and supersaturation PDFs and Lagrangian statistics?

Section 2 theoretically discusses the contribution of phase change to water vapor and temperature variances, their covariance, and supersaturation variability. It also introduces a subgrid model for these contributions, extending the subgrid scheme for without cloud conditions proposed in Chandrakar et al. (2022). Section 3 discusses the DNS model and simulation setup. The results section discusses the Lagrangian (sec. 4a) and Eulerian (4b) statistics of water vapor, temperature, and supersaturation and evaluates the model introduced in sec. 4d. Finally, sec. 5 and 6 discuss the implications of our findings for subgrid supersaturation modeling and the primary conclusions of this study.

Subgrid-scale turbulent variability at different grid cell locations in an Eulerian domain is needed to represent the subgrid microphysical interactions in Eulerian or Lagrangian microphysics models. As discussed in the introduction, such Eulerian subgrid information can be used to drive the magnitude of the Weiner forcing term in Lagrangian stochastic models of scalar fluctuations.

Below, an Eulerian subgrid supersaturation model is discussed that includes the contribution from the phase change in clouds, extending the model presented in Chandrakar et al. (2022).

The variance of subgrid supersaturation (𝑠) fluctuations can be expressed as a function of the variances and covariances of potential temperature (𝜃), water vapor mixing ratio (𝑞 𝑣 ), and pressure (𝑝) or vertical velocity (𝑤) (Kulmala et al. 1997;Siebert and Shaw 2017;Chandrakar et al. 2020Chandrakar et al. , 2022)):

where denotes the spatial filtering operation for different variables, 𝛽 = 𝐿/(𝑅 𝑣 T), T is the resolved scale temperature, 𝑅 𝑣 and 𝑅 𝑑 are the gas constants for water vapor and air, 𝐿 is the latent heat of condensation, 𝐶 𝑝 is the specific heat capacity at constant pressure, 𝑔 is the gravitational acceleration, 𝑆 = q𝑣 /𝑞 𝑣𝑠 ( θ), 𝑞 𝑣𝑠 is the saturation vapor mixing ratio, and 𝜏 𝑡 is the eddy turnover time or the turbulent mixing timescale. The approximate sign here is due to neglecting terms higher than the second order. The budget equations for subgrid water vapor and potential temperature variances and covariances are derived elsewhere (Deardorff 1974a;Wyngaard et al. 1978;Chandrakar et al. 2022) but are presented in Appendix A for convenience. These budget equations can be solved to obtain water vapor and temperature variances and their covariance, determining the supersaturation variance. This model is evaluated against DNS in moist convection without cloud formation in Chandrakar et al. (2022).

In the subgrid model and simulations presented in Chandrakar et al. (2022), the diabatic contributions of condensation/evaporation (the last terms in Eqs. A1, A2, and A3) to the prognostic scalar variances and covariance were neglected because there were no droplets and hence no water phase change. However, in cases with droplets, the diabatic terms relax the saturation ratio toward unity. In our supersaturation fluctuation model, the impact of condensation can be introduced as an additional sink term (-(𝑠s)/𝜏 𝑐 ), and it reduces the supersaturation variance at the rate of -2( 𝑠𝑠 -s2 )/𝜏 𝑐 (assuming a homogeneous distribution of cloud droplets). Here, 𝜏 𝑐 is the phase relaxation timescale of cloud droplets and depends on the first moment of the droplet number concentration distribution. Past studies (e.g., Chandrakar et al. 2016Chandrakar et al. , 2017;;MacMillan et al. 2022) also suggested that increased aerosol concentration (and a corresponding increase in the phase relaxation timescale) suppresses supersaturation variability. In the budget equations for subgrid water vapor and temperature variances and covariance, the contribution of phase change depends on 𝑠𝑞 𝑣 , 𝑠𝜃, and 𝜏 𝑐 . These diabatic terms are represented as:

where 𝑞 𝑙 is the condensation/evaporation rate, Π = (𝑃/𝑃 𝑜 ) 𝑅 𝑑 /𝐶 𝑝 is the Exner function, 𝑃 𝑜 is a reference surface air pressure, 𝑎 1 ( θ) = 𝑞 -1 𝑣𝑠 (𝜕𝑞 𝑣𝑠 /𝜕𝜃) is a function of resolved scale θ, p, and some thermodynamics constants. 𝜏 * 𝑐 = 𝜏 𝑐 /𝜁 = (4𝜋𝜌 -1 𝑎 𝜌 ℓ 𝜉 𝐼) -1 is the phase relaxation timescale divided by an additional thermodynamic factor 𝜁 = (𝑞 𝑣𝑠 + 𝐿 2 /(𝑅 𝑣 𝑇 2 𝐶 𝑝 )) -1 to separate the effect of change in air temperature from the latent heat release in the phase relaxation timescale (separately accounted through the heat equation). 𝜌 𝑎 and 𝜌 ℓ are the air and liquid water density, respectively. 𝜉 ( T, q𝑣 ) is the thermodynamic growth factor in the droplet growth equation and depends on the mean temperature and water vapor (see eq. 7.17 in Rogers and Yau 1996). 𝐼 = ∫ 𝑅𝑛 𝑑 (𝑅)𝑑𝑅 is the first moment of the droplet size distribution (integral radius), where 𝑛 𝑑 (𝑅) is the droplet number concentration distribution and 𝑅 is the droplet radius. Note here all diabatic contributions are simplified as a function of water vapor and temperature variances and their covariance. Thus, the budget equations presented in Appendix A are fully closed.

The above expressions are obtained using the droplet condensation growth equation ( 𝑞 𝑙 ≈ 4𝜋𝜌 ℓ 𝜌 -1 𝑎 𝜉 𝐼 𝑠) and a Taylor series expansion of saturation ratio fluctuations, including first-order terms only (similar to the one used for Eq. 2). Note these expressions neglect the terms involving covariance of the integral radius with supersaturation, water vapor mixing ratio, and temperature. These terms might be important near cloud boundaries where entrainment could lead to significant inhomogeneity due to in-cloud activation, evaporation, and dilution. For simplicity, we neglect these contributions here. Future work will investigate the impact of these terms in the context of nonuniform mixing. Subgrid condensation introduces a positive heat flux and negative water vapor flux at small scales, with the opposite sign of fluxes for subgrid evaporation. Thus, the diabatic terms decrease the covariance of water vapor and temperature fluctuations but increase the fluctuation magnitude for both scalars individually at small scales. The phase change contribution of the simplified subgrid model proposed here, Eqs. A1, A2, and A3, is evaluated in section 4d.

The moist Rayleigh-Bénard DNS setup follows earlier work presented in Chandrakar et al. (2022) using Cloud Model 1 (CM1, Bryan and Fritsch 2002). The setup is similar to the Pi Cloud Chamber (Chang et al. 2016;Chandrakar et al. 2016), with domain size = 1.08 m × 1.08 m × 1 m, Rayleigh number = 1.52 × 10 9 , Prandtl number = 0.72, and Schmidt number = 0.62, except the aspect ratio (width to height ratio), which is unity here instead of two for the chamber. The Taylor microscale Reynolds number (𝑅𝑒 𝜆 ) in the convection core for these simulations is around 96-123.

CM1 solves the compressible Boussinesq equation set (Boussinesq approximation for buoyancy in the momentum equation with a prognostic pressure equation) using a three-step Runge-Kutta time integration method with a 5th order advection scheme. Scalar advection uses a weighted essentially nonoscillatory (WENO) scheme. An adaptive timestep with maximum CFL limit of 0.9 is used for integrating the governing equation set. The compressible solver employed here uses Klemp-Wilhelmson time-split steps for acoustic terms. A no-slip velocity boundary condition is used for all walls. The top and bottom boundaries use a fixed water vapor and temperature boundary condition (saturated at 282 K and 294 K, respectively). The sidewall conditions vary for different cases to mimic relative scalar fluxes during entrainment and to match realistic conditions in the Pi Cloud Chamber. Note that the purpose of the current simulations is not to exactly mimic the entrainment water vapor and heat fluxes for a given turbulent intensity in atmospheric clouds, but rather to understand the impact of relative fluxes of water vapor and temperature on the scalar variances and covariance and, ultimately, supersaturation variability. Table 1 summarizes different simulation cases characterizing different degrees of entrainment from the side in moist and cloudy conditions. Please refer to Chandrakar et al. (2022) for more details about the model and the simulation setup.

For the current study, we extended the moist (with water vapor, but without aerosol or cloud particles) simulations from Chandrakar et al. (2022) after introducing NaCl aerosol particles of a fixed size (130 nm) and allowing the system to reach a steady-state condition. NaCl aerosols are also continuously injected at a constant rate (33 cm -3 min -1 ) every 5 s interval to allow a balance between particle loss at boundaries and injection. The injected aerosol size and rate are similar to a recent experimental study in the Pi Cloud Chamber with one of the intermediate aerosol injection rates (see Table 1 in Prabhakaran et al. 2020). At the start of the simulations, 18.7 million Lagrangian computational particles (with an extra 110% buffer to allow injection) are used. For the current low aerosol condition and without droplet collision-coalescence, this number is expected to be sufficient for representing cloud microphysical interactions and droplet statistics.

A Lagrangian microphysics scheme based on Shima et al. (2009) (with some modification) is used here. Droplet condensation and evaporation are calculated from the grid-scale water vapor and temperature fields, and tendencies for water vapor and temperature at each grid cell are passed to the Eulerian dynamical core after each condensation/evaporation step. All particles are assumed to have no thermal inertia since cloud droplets equilibrate with their neighboring environment significantly faster (O (10 -5 s) ) than the timescale of smallest eddies. The liquid water content at each grid is also provided to the dynamical core for its buoyancy calculation. Each Lagrangian computational point particle carries information on solute mass, and the scheme explicitly solves the droplet condensation/evaporation equation with curvature and solute terms. Köhler theory is used for droplet activation. For particle transport, the following transport equation is used assuming the Stokes limit:

where ì 𝑑𝑣 𝑝 is the particle velocity vector, ì 𝑢 is the fluid velocity linearly interpolated to the location of a particle, ì 𝑔 is the gravitational acceleration vector, and 𝜏 𝑑 = 2𝑅 2 𝜌 ℓ /(9𝜇) is the droplet response time (𝑅 is the droplet radius, and 𝜇 is the dynamic viscosity of air). Both the droplet condensation/evaporation and transport equations are solved using an implicit method with the same timestep as the dynamical core.

These simulations also include an additional five thousand passive point particles to compare Lagrangian statistics of passive particles with inertial droplets. The Lagrangian data of two to five thousand randomly positioned droplets and passive particles is output for the last 2 min of the simulation at 0.1 s intervals (less than half of the Kolmogorov timescale, 𝜏 𝜂 ) for calculating turbulence statistics. The three-dimensional Eulerian fields are output every 30 s (order of the large-scale circulation time (Niedermeier et al. 2018)), and the Eulerian statistics are output every 20 s. The turbulent statistics presented here are calculated over a period of several hundred times the free-fall buoyancy timescale:

where 𝐻 is the height of the domain, 𝑔 is the gravitational acceleration, 𝛼 is the thermal expansion coefficient, and Δ𝑇 𝑣 is the virtual temperature difference between the top and bottom boundaries.

Table 1. Summary of moist and cloudy DNS cases presented in this article. The Rayleigh number for these cases is 1.52 × 10 9 . The isotropic grid spacing (Δ) is normalized based on the grid resolution criteria for the Rayleigh-Bénard convection in Grötzbach (1983). Here, ⟨𝜂⟩ is the average Kolmogorov length scale in the convection core from DNS runs, and 𝑞 𝑣 𝑠 (𝑇 𝑠 ) is the saturation vapor mixing ratio at the sidewall temperature.

Note -⟨ * ⟩ denotes ensemble averaging.

# Case Number of grid points Grid resolution Simulation type Sidewall condition 1. DNS-AD 540 × 540 × 500 𝜋 ⟨ 𝜂⟩ / Δ = 2.9 DNS (no cloud) Adiabatic sidewalls 2. DNS-S100 540 × 540 × 500 𝜋 ⟨ 𝜂⟩ / Δ = 2.9 DNS (no cloud) mean temperature and 100 % of 𝑞 𝑣𝑠 (𝑇 𝑠 ) 3. DNS-S95 540 × 540 × 500 𝜋 ⟨ 𝜂⟩ / Δ = 2.9 DNS (no cloud) mean temperature and 95 % of 𝑞 𝑣𝑠 (𝑇 𝑠 ) 4. DNS-AD-CLD 540 × 540 × 500 𝜋 ⟨ 𝜂⟩ / Δ = 2.9 DNS (cloud) Adiabatic sidewalls 5. DNS-S100-CLD 540 × 540 × 500 𝜋 ⟨ 𝜂⟩ / Δ = 2.9 DNS (cloud) mean temperature and 100 % of 𝑞 𝑣𝑠 (𝑇 𝑠 ) 6. DNS-S95-CLD 540 × 540 × 500 𝜋 ⟨ 𝜂⟩ / Δ = 2.9 DNS (cloud) mean temperature and 95 % of 𝑞 𝑣𝑠 (𝑇 𝑠 ) 13 Accepted for publication in Journal of the Atmospheric Sciences. DOI 10.1175/JAS-D-22-0256.1. Brought to you by MICHIGAN TECHNOLOGICAL UNIVERSITY | Unauthenticated | Downloaded 08/07/23 05:29 PM UTC There is also variability in the lifetime of particles depending on their location and growth history, consistent with previous DNS of the Pi Chamber (MacMillan et al. 2022). 20 10 0 10 1 10 2 -1 10 1 dq 16 Accepted for publication in Journal of the Atmospheric Sciences. DOI 10.1175/JAS-D-22-0256.1. Brought to you by MICHIGAN TECHNOLOGICAL UNIVERSITY | Unauthenticated | Downloaded 08/07/23 05:29 PM UTC -0.05 0 0.05 s' -4 -3 -2 -1 PDF -0.05 0 0.05 s' -4 -3 -2 -1 PDF -0.05 0 0.05 s' -4 -3 -2 -1 PDF -0.05 0 0.05 s' 10 -4 10 -3 10 -2 10 -1 PDF -0.05 0 0.05 s' 10 -4 10 -3 10 -2 10 -1 PDF -0.05 0 0.05 s' 10 -4 10 -3 10 -2 10 -1 PDF (a) DNS-AD s =0.4 % Skewness = -7.4 Kurtosis = 80.3 (b) DNS-AD s = 0.33% Skewness =3.1 Kurtosis = 16.6 (c) DNS-S100 s =0.56% Skewness = -0.46 Kurtosis = 6.1 (d) DNS-S100 s = 0.32% Skewness = 3.77 Kurtosis = 22.2 (f) DNS-S95 s = 1.08% Skewness = -0.96 Kurtosis = 5.2 (e) DNS-S95 s = 1.34% Skewness = -0.87 Kurtosis = 4.62 Before Cloud After Cloud Droplets R > 5 m: s = 0.22% Skewness = 3.3 Kurtosis = 20.3 Fig. 4. Lagrangian PDFs of supersaturation before (left) and after (right) cloud formation. Red lines show Gaussian fits. 17 Accepted for publication in Journal of the Atmospheric Sciences. DOI 10.1175/JAS-D-22-0256.1. Brought to you by MICHIGAN TECHNOLOGICAL UNIVERSITY | Unauthenticated | Downloaded 08/07/23 05:29 PM UTC -6 -4 -2 0 2 4 6 ds( ) ds 2 -1/2 10 -4 10 -3 10 -2 10 -1 PDF / = 0.44 Gaussian Fit / = 0.89 / = 1.78 / = 4.45 / = 8.90 -6 -4 -2 0 2 4 6 ds( ) ds 2 -1/2 10 -4 10 -3 10 -2 10 -1 PDF / = 0.42 Gaussian Fit / = 0.84 / = 1.68 / = 4.19 / = 8.38 -6 -4 -2 0 2 4 6 ds( ) ds 2 -1/2 10 -4 10 -3 10 -2 10 -1 PDF / = 0.39 Gaussian Fit / = 0.78 / = 1.55 / = 3.89 / = 7.77 -6 -4 -2 0 2 4 6 ds( ) ds 2 -1/2 10 -4 10 -3 10 -2 10 -1 PDF / = 0.47 Gaussian Fit / = 0.93 / = 1.87 / = 4.67 / = 9.35 -6 -4 -2 0 2 4 6 ds( ) ds 2 -1/2 10 -4 10 -3 10 -2 10 -1 PDF / = 0.42 Gaussian Fit / = 0.84 / = 1.67 / = 4.18 / = 8.36 -6 -4 -2 0 2 4 6 ds( ) ds 2 -1/2 10 -4 10 -3 10 -2 10 -1 PDF / = 0.42 Gaussian Fit / = 0.85 / = 1.69 / = 4.23 / = 8.47 After Cloud Before Cloud (a) DNS-AD (b) DNS-AD (c) DNS-S100 (f) DNS-S95 (d) DNS-S100 (e) DNS-S95 Fig. 5. PDFs of Lagrangian supersaturation difference at different lags (𝜏/𝜏 𝜂 ) before (left) and after (right) cloud formation.

Lagrangian statistics of turbulent velocity fields have been extensively studied using measurements and DNS, especially in the context of developing Lagrangian stochastic models of turbulent transport (e.g., Pope 1994;Yeung 2002;Biferale et al. 2006Biferale et al. , 2008)). Although Lagrangian scalar statistics were less of a focus in past investigations, they are critical for understanding cloudturbulence interactions (Paoli and Shariff 2009). Supersaturation, which is a nonlinear combination of two scalars (temperature and water vapor), brings another level of complication in investigating and modeling turbulence from a Lagrangian perspective. In this section, we discuss different Lagrangian statistics that provide the basis for subgrid modeling of supersaturation variability and its interactions with cloud particles. Note that the turbulent Rayleigh-Bénard convection is inherently non-isotropic (except near the convection core). However, most of the scaling laws of turbulent statistics (e.g., the Kolmogorov-Obukhov-Corrsin scaling) assume isotropic turbulence. Thus, care should be taken while interpreting the presented results.

Figure 2 shows Lagrangian autocorrelation functions of water vapor, temperature, and supersaturation before and after cloud formation for different sidewall boundary conditions. The integral or autocorrelation timescale is calculated by integrating the autocorrelation function from zero lag to the first zero crossing. These autocorrelation functions are close to exponential, which is expected for velocity and passive scalars in turbulent flow (Tennekes 1979). The supersaturation, a joint scalar, also displays a nearly exponential Lagrangian autocorrelation, thus justifying a key assumption of Lagrangian stochastic models (since this is one of the criteria for a Markovian process). Except for the driest sidewall condition, the water vapor and temperature Lagrangian integral timescales are nearly the same. For the driest sidewall condition, the water vapor integral timescale is somewhat longer without cloud compared to the other cases, which shows the impact of a significant sidewall vapor flux. However, when cloud is present the water vapor integral timescale is nearly the same as that for the other sidewall conditions.

The most significant impact of sidewall fluxes appears in the supersaturation autocorrelation.

The supersaturation autocorrelation timescale increases with an increase in sidewall fluxes in the absence of cloud. It is worth noting that the supersaturation autocorrelation timescale differs from the Lagrangian timescales for water vapor and temperature, and is lower for adiabatic sidewalls and higher for non-adiabatic sidewalls before cloud formation. Thus, assuming the supersaturation autocorrelation timescale is the same as individual scalars or velocity may not be reasonable, especially in the presence of a relative external flux.

The presence of cloud suppresses the supersaturation correlation time for non-adiabatic sidewalls.

For adiabatic sidewall conditions, cloud formation slightly increases the correlation time compared to that before cloud. Interestingly, this makes the correlation timescale for the adiabatic case larger than for the saturated sidewall condition. In the adiabatic case before cloud formation, water vapor and temperature fields are highly correlated due to the absence of significant relative fluxes. When cloud forms, relative fluxes from condensation decreases the correlation between the two, increasing the supersaturation correlation timescale. However, condensation/evaporation acts as a high-pass filter of supersaturation fluctuations, suppressing extreme decorrelation when temperature and water vapor are already highly decorrelated. Lagrangian supersaturation autocorrelation functions for cloud droplets with radius 𝑅 > 5 𝜇m are slightly suppressed compared to the passive (massless) particle trajectories. For these cloud droplets, gravitational settling acts as a low-pass filter and reduces the correlation of droplets with supersaturation field slightly. This impact on correlation will likely be further enhanced for larger cloud and rain drops. However, for this size range, drop collision-coalescence might be more important than condensation for drop growth, and the suppression of supersaturation timescale may not be critical. The autocorrelation functions of the Lagrangian temporal gradients of water vapor, temperature, and supersaturation (not shown) are extremely narrow (zero crossing at < 1 s), thus justifying the assumption of delta correlated Wiener forcing in Lagrangian stochastic models of these quantities. A similar scaling is used to normalize the Lagrangian supersaturation structure functions as that for independent scalars (𝜏) but with a different scaling coefficient (2𝜎 2 𝑠 /𝜏 𝑡 𝑠 , where 𝜏 𝑡 𝑠 ) is the Lagrangian supersaturation autocorrelation timescale). This scaling coefficient can be obtained from the exponential autocorrelation function in a steady-state condition, and is used for simplicity since the dissipation rate is not defined for supersaturation. An equivalent parameter could be derived for supersaturation but would depend on complex high-order covariances of scalar gradients. As shown in Fig. 3(e-f), supersaturation displays inertial range scaling similar to the independent scalars. However, the supersaturation second-order structure functions show a slightly wider plateau (or peak) region than those for the independent scalars, extending toward smaller and larger lags. The positive slope in the dissipation range is also weaker. Gotoh et al. (2021) also showed an Eulerian inertial range supersaturation spectrum (𝑘 -5/3 ) consistent with our second-order Lagrangian structure functions, but they treated supersaturation similar to independent scalars. The Eulerian one dimensional energy, independent scalar, and supersaturation spectra away from boundaries (not shown) also display a power-law region with a slope close to 𝑘 -5/3 scaling for all cases presented here.

Figure 3 also compares the structure functions for different sidewall conditions. Sidewall fluxes do not have a significant impact on the shape of structure functions for water vapor and temperature. This is expected since additional fluxes should only affect the scalar variability magnitude but not how the scalars mix by turbulence. However, this is not the case for supersaturation because of its nonlinear dependence on the independent scalars. In the absence of cloud, the peak of the supersaturation structure functions shifts towards a larger value with increasing sidewall relative fluxes. Moreover, the negative slope becomes shallower at larger lags (larger 𝜏 in the energy containing range), particularly for the case with subsaturated sidewalls (DNS-S95). For this case, the positive slope at small lag is also steeper than in the other two runs. Since the sample trajectories do not exclude the regions near boundaries, some of these effects at large lags could be associated with the influence of boundary fluxes.

With cloud, condensation and evaporation do not significantly affect the shape of the water vapor and temperature structure functions relative to non-cloudy conditions. In the current cases, the droplet population introduces fluxes of water vapor and latent heat from phase change at a timescale (the phase relaxation timescale 𝜏 𝑐 /𝜏 𝜂 = 37, 30, and 140 for the adiabatic, S100, and S95 sidewall conditions) much longer than the timescale of inertial range eddies. Thus, the impact of phase change on the water vapor and temperature fluctuations is expected to be minimal. However, supersaturation is more sensitive to the relative fluxes of temperature and water vapor from phase change (not just their individual magnitudes) owing to its nonlinear dependence. Consequently, condensation and evaporation suppress the far right tail of the supersaturation structure function for the saturated sidewall case towards a stronger slope. But for the adiabatic sidewalls, it enhances the relative magnitude of the right tail (i.e., weaker slope). For the subsaturated sidewall condition, the right tail is not significantly suppressed (slope is not increased) since the phase relaxation timescale is substantially longer.

Droplet inertia could also impact the scalar and supersaturation structure functions, but most cloud droplets have an inertial response timescale much shorter than the Kolmogorov timescale (i.e., 𝑆𝑡 ≪ 0.1). Thus, no significant impact is expected from cloud droplet inertia. However, the impact of gravitational settling cannot be ignored. The ratios of the gravitational settling velocity (𝑣 𝑠𝑡 ) and the vertical flow velocity scale ((⟨𝑤 ′ 2 ⟩) 1/2 ) in the current case for 5 𝜇m, 10 𝜇m, and 20 𝜇m radius droplets are 0.05, 0.19, and 0.76, respectively. Thus, droplet settling causes slight decorrelation with flow trajectories. The water vapor and temperature structure functions for drops with 𝑅 > 5 𝜇m have substantially lower magnitudes at all scales compared to those for passive parcels, as shown in Fig. 3b,d, due to droplet settling. However, the shape remains the same and thus there is little difference in the supersaturation structure functions for drops compared to passive parcels (Fig. 3f).

The shape of the Lagrangian supersaturation PDF in turbulent convection is another important factor for understanding and parameterizing the effect of turbulent fluctuations on droplet activation and condensation growth. Figure 4 shows Lagrangian supersaturation fluctuation PDFs averaged over all passive particles before and after cloud formation. These PDFs of Lagrangian fluctuations (𝑠 ′ = 𝑠(𝑡) -⟨𝑠(𝑡)⟩ 𝑡 ) are calculated for each sampled particle over 1-1.5 min of the sampling time and averaged over all particles. Here, ⟨ * ⟩ 𝑡 denotes averaging along a Lagrangian trajectory. Note that these PDFs of randomly sampled trajectories also include trajectories passing through the regions near boundaries, thus sampling a few extreme values of scalars from the viscus layers. Separating those trajectories and maintaining a long enough sampling time for a sufficient number of sampled particles was challenging. Nevertheless, the average time spent by randomly sampled particles near boundaries is expected to be significantly lower compared to the main region away from the boundaries. For adiabatic sidewalls, the supersaturation is highly negatively skewed before cloud formation, although the water vapor and temperature PDFs (not shown) are only slightly positively skewed. Moreover, the kurtosis value of ∼80 for adiabatic sidewalls is significantly higher than that for a Gaussian PDF (kurtosis = 3) or water vapor and temperature individually (21-23, not shown).

For the cases with non-adiabatic sidewalls, the supersaturation PDFs have greatly reduced negative skewness and kurtosis before cloud formation. The value of kurtosis (4.6-6.1) is also much closer to Gaussian than the water vapor or temperature PDFs (25-28, not shown), and it decreases with an increase in sidewall relative fluxes.

Cloud formation has a dramatic impact on the Lagrangian supersaturation PDFs and changes the sign of skewness for the adiabatic and saturated sidewall cases. It is positive for these cases in contrast to the negative skewness before cloud. The kurtosis decreases for the adiabatic case but increases for the saturated sidewalls. The case with subsaturated sidewalls has only minor differences in skewness and kurtosis before and after cloud formation. For all cases, the magnitude of supersaturation fluctuations decreases somewhat with cloud formation. Moreover, along the trajectories of large cloud droplets (𝑅 > 5 𝜇m), the supersaturation variability is suppressed further.

Note that the Lagrangian statistics of large cloud droplets are only shown for the DNS-AD-CLD case since the droplets are larger in this case. Thus, the inertial and settling impacts are expected to be more significant than the other two cases. Moreover, the probability of sampling the large cloud droplets with a sufficiently long sampling time is also significantly low for the other two cases.

The PDF shape of the Lagrangian difference (or increment) is another critical foundation for the stochastic models (typically assumed Gaussian). The Lagrangian difference is the difference in a fluid quantity along particle trajectories at different time lags: 𝑑𝑠(𝜏) = 𝑠(𝑡 + 𝜏) -𝑠(𝑡). These PDFs are calculated similarly to the Lagrangian fluctuation PDFs for each sampled particle over 2 min of the sampling time and averaged over all particles. PDFs of Lagrangian fluid acceleration and Lagrangian velocity differences at different lags in turbulent flow has been extensively studied in the past to support the formulation of Lagrangian stochastic models of particle transport (Pope 1994;La Porta et al. 2001;Yeung 2002;Reynolds et al. 2005;Biferale et al. 2006;Brown et al. 2009). Figure 5 shows average PDFs (from the passive particles) of Lagrangian supersaturation difference at various lags ranging from sub-Kolmogorov to inertial range (𝜏/𝜏 𝜂 = 0.39 -9.35).

The shape of these PDFs for different lags is nearly the same. They are close to Gaussian near the mode, but the tails deviate significantly from Gaussian (closer to an exponential shape).

This shape of Lagrangian supersaturation difference PDFs (analogous to the shape of Lagrangian supersaturation gradient PDFs at small lags), with significant intermittency, is reminiscent of the PDF shape of Lagrangian fluid acceleration for particle transport observed in fully developed turbulence (La Porta et al. 2001;Brown et al. 2009). Note that Lagrangian supersaturation gradient PDFs are analogous to acceleration PDFs but in the space of 𝑅 2 for droplet growth instead of fluid velocity. The PDF shape of the Lagrangian difference for independent scalars (not shown) is similar to the supersaturation difference PDFs presented here. PDFs of Lagrangian differences are more symmetric than the supersaturation PDFs themselves. However, with cloud formation, the tails become slightly asymmetric (positively skewed). This is most apparent in the adiabatic sidewall case. Eulerian scalar statistics provide a different perspective than the Lagrangian statistics and are also crucial for understanding and modeling the interactions between cloud droplets and scalar fields. This section discusses some important Eulerian statistics of the scalars and supersaturation before and after cloud formation for the different sidewall boundary conditions. Figure 6 shows horizontally-averaged profiles (excluding 40 grid points along all four sidewalls) of the standard deviation of water vapor, temperature, supersaturation, and water vapor-temperature covariance. As shown in Chandrakar et al. (2022), there is no significant difference in water vapor and temperature standard deviations among the simulations in the convection core (the region away from sidewalls).

With cloud formation, the changes in water vapor and temperature variability in the convection core are minor. However, the covariance of water vapor and temperature is sensitive to both sidewall relative fluxes and cloud formation. Chandrakar et al. (2022) showed the response of the covariance magnitude to relative sidewall fluxes. Here, we show that cloud formation increases the normalized covariance in the presence of relative sidewall fluxes. However, it decreases the cross-correlation of water vapor and temperature for the adiabatic sidewall case, where these scalars are highly correlated in the absence of clouds. This impact on covariance is also reflected in the standard deviation of supersaturation that depends on the covariance magnitude. Supersaturation standard deviation decreases for the cases with relative sidewall fluxes but increases for the adiabatic sidewalls after cloud formation, consistent with the response of the normalized covariance.

The Eulerian PDF of supersaturation is also important for modeling condensation/evaporation and feedback to the dynamics. It is also interesting from the point of fundamental fluid dynamics; how does a nonlinear transformation of scalar fields (supersaturation in the current case) in turbulent flow affect the Eulerian PDF of the transformed scalar? In homogeneous and isotropic turbulence, scalar PDFs are expected to be close to Gaussian with slight deviation near the tails representing intermittency. However, the presence of a large-scale gradient in a scalar field causes deviation from the Gaussian PDF shape (Shraiman and Siggia 2000). In Rayleigh-Bénard convection, the temperature gradient between the top and bottom boundaries introduces such a large-scale scalar gradient. Figure 7 shows Eulerian PDFs of temperature fluctuations in the convection core (50 or more grid points away from all boundaries) for different sidewall conditions. The Gaussian reference curve is obtained from fitting the PDFs data to a Gaussian function. In all cases, the temperature PDFs significantly deviate from Gaussian shape. However, they are nearly symmetric (close to zero skewness). The PDF shape is close to an exponential with kurtosis of 11-13 (compared to 3 for a Gaussian PDF). The presence of the cloud does not significantly impact the overall PDF shape. Previous studies also reported nearly exponential PDF tails of temperature fluctuations in the core of Rayleigh-Bénard convection (Castaing et al. 1989;Gollub et al. 1991;Wunsch and Kerstein 2005) as well as Rayleigh-Bénard convection with water vapor (Chandrakar et al. 2020).

The reason for the deviation of the tails from Gaussian is the injection of fluid from the top and bottom boundary layers (Castaing et al. 1989;Wunsch and Kerstein 2005). It is also important to note that the sidewall boundary condition has little influence on the temperature PDF shape since the net sidewall heat flux is relatively small. A similar situation may occur in atmospheric clouds where entrainment introduces large-scale gradients in scalar fields that are subsequently smoothed out due to mixing. Thus, a non-Gaussian shape of scalar PDFs could also be present in clouds.

The Eulerian water vapor fluctuation PDFs are shown in Fig. 8. For adiabatic sidewalls, the water vapor PDF shape is nearly the same as the temperature PDF (nearly zero skewness and exponential tails), similar to a past study of Rayleigh-Bénard convection with water vapor (Chandrakar et al. 2020). However, for non-adiabatic sidewalls, the water vapor PDFs become more positively skewed. The skewness magnitude is higher for DNS-S95 (and 90% saturated sidewalls, not shown) compared to DNS-S100. This positive skewness is due to the increased water vapor flux at the lower boundary for non-adiabatic sidewalls (see Chandrakar et al. 2022).

It leads to a higher concentration of water vapor in the plumes ejected from the bottom boundary layer, enhancing positive fluctuations in the convection core. Tails of these PDFs have a similar deviation from Gaussian as the temperature PDFs, and their kurtosis is between 10-13. After cloud formation, the PDF shape remains nearly the same, but the standard deviation slightly decreases.

The standard deviation of water vapor fluctuations is somewhat higher for DNS-S95-CLD than for DNS-AD-CLD and DNS-S100-CLD. The skewness value is also slightly higher with cloud for the adiabatic case but remains similar for the other cases compared to before cloud. Accepted for publication in Journal of the Atmospheric Sciences. DOI 10.1175/JAS-D-22-0256.1. Brought to you by MICHIGAN TECHNOLOGICAL UNIVERSITY | Unauthenticated | Downloaded 08/07/23 05:29 PM UTC 1 2 3 Case # -1 -0.5 0 0.5 1 Contribution to s' 3 10 -5 I II III IV 1 2 3 Case # -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 Contribution to s' 3 s -3 I II III IV 1 2 3 Case # -1 -0.5 0 0.5 1 Contribution to s' 3 10 -5 I (cloud) II (cloud) III (cloud) IV (cloud) 1 2 3 Case # -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 Contribution to s' 3 s -3 I (cloud) II (cloud) III (cloud) IV (cloud) (a) (b) (c) (d) Fig. 10. Different terms marked (I -IV) in Eq. 7 contributing to supersaturation skewness. The left plots (a and c) show contributions to the third moment of supersaturation fluctuations as in Eq. 7, and the right plots (b and d) show the normalized contributions to supersaturation skewness. The case numbers are the same as in Table 1 (Case 1: DNS-AD, Case 2: DNS-S100, Case 3: DNS-S95). Results for non-cloud conditions are shown in the top plots (a and b), and those corresponding to cloudy conditions are presented at the bottom plots (c and the DNS-AD case; these PDFs are much closer to a Gaussian distribution. This is a key finding because it supports the assumption of a Gaussian PDF in the subgrid condensation models in more realistic cases where relative water vapor and temperature fluxes are usually present.

After cloud formation, the Eulerian supersaturation PDFs become positively skewed for the adiabatic and saturated sidewall cases, consistent with the Lagrangian PDFs. However, they becomes slightly negatively skewed for the subsaturated sidewall case. The kurtosis value is also higher with cloud for the saturated and subsaturated sidewall cases compared to before cloud.

Moreover, the supersaturation standard deviation decreases after cloud formation for non-adiabatic sidewalls and increases for the adiabatic case. This is also seen in the vertical profiles of 𝜎 𝑠 in Fig. 6d.

To investigate the supersaturation PDF skewness further, we derived an expression for the supersaturation skewness in terms of individual scalar fluctuations. It is a function of the mean temperature and water vapor mixing ratio, skewness of temperature and water vapor fluctuations, ⟨𝑇 ′2 𝑞 ′ 𝑣 ⟩, ⟨𝑇 ′ 𝑞 ′2 𝑣 ⟩, and higher-order fluctuation terms (that are neglected):

Here, ⟨ * ⟩ denotes ensemble averaging. This expression for supersaturation skewness (without normalization by 𝜎 3 𝑠 ) is obtained from a Taylor series expansion and neglecting the higher-order contributions (see Appendix B for details). Contributions from the skewness of individual scalars are opposite of one another; the water vapor skewness term is positive, and the temperature skewness term is negative. The third and fourth covariance terms also have opposite contributions to the supersaturation skewness. The last covariance term has a negative contribution.

Figure 10 shows the magnitude of different terms from Eq. 7. For DNS-AD, the absolute values of all the terms are small, and they mostly cancel each other out. However, the normalized contributions (normalized by 𝜎 3 𝑠 ) are the largest among all cases since 𝜎 𝑠 is very small. Thus, the assumption of negligible higher-order terms in the derivation of Eq. 7 is invalid for DNS-AD. A significant negative skewness of the supersaturation PDF for this case is caused by these neglected higher-order terms. For more realistic cases involving sidewall fluxes (before cloud), the absolute magnitude of different terms increases with decreasing sidewall saturation level (i.e., more negative flux of water vapor). The most significant contribution in the cases with sidewall fluxes is the term associated with ⟨𝑞 ′2

𝑣 𝑇 ′ ⟩. The water vapor skewness term is positive and is the secondlargest contributor. The term with ⟨𝑞 ′ 𝑣 𝑇 ′2 ⟩ is similar in magnitude to the water vapor skewness term. The temperature skewness term is almost zero. In cases with significant external fluxes, the ⟨𝑞 ′ 𝑣 𝑇 ′2 ⟩ contribution and the temperature skewness terms should be more significant. However, similar to the current simulations, the normalized values of different skewness terms depend on the supersaturation fluctuation magnitude. Moreover, the relative magnitude and sign of external heat and moisture fluxes in the system will determine whether the net supersaturation skewness will be zero, positive or negative. Thus, entrainment in the atmospheric boundary layer, which tends to produce negatively skewed water vapor and positively skewed temperature fluctuations (Deardorff 1974b;Mahrt 1991;Mellado et al. 2017), is expected to increase skewness in the supersaturation PDF (considering only the scalar skewness terms). The other third-order correlation terms (the third and fourth terms) would also be affected by cloud-top entrainment. In the case of lateral entrainment in cumulus clouds, entrained environmental air is typically colder and dryer than the convection core (Katzwinkel et al. 2014). Thus, lateral mixing of environmental air into a cloud is expected to produce more negative fluctuations and skewness in water vapor and temperature PDFs.

Hence, they would be expected to compensate one another to produce a less skewed supersaturation PDF.

In the presence of cloud, the skewness of water vapor fluctuations and other third-order covariances (⟨𝑞 ′2 𝑣 𝑇 ′ ⟩ and ⟨𝑞 ′ 𝑣 𝑇 ′2 ⟩) remain more significant contributors to supersaturation skewness than the temperature skewness. However, a slight negative skewness of temperature fluctuations (a positive second term) after cloud formation seems to be important for explaining the change in sign of the supersaturation skewness from negative to positive for the adiabatic case (keeping in mind that the neglected higher order terms may also be important for this case). For the saturated sidewall case, a small decrease in the magnitude of the fourth term (⟨𝑇 ′ 𝑞 ′2

𝑣 ⟩) appears to explain the change in supersaturation skewness from nearly zero to positive after cloud formation. For the subsaturated sidewall case, all four terms undergo small changes with cloud formation such that the supersaturation skewness changes from near zero to slightly negative. Subgrid modeling of water vapor and temperature variances, their covariance, and supersaturation variance in cloudy conditions requires modeling the contributions from phase change (condensation and evaporation). Section 2b introduced a simplified representation (Eqs. 4, 3, 5) of these contributions in the budget equations for water vapor and temperature variance and covariance. An a priori test is used to evaluate the proposed formulations. We use water vapor and temperature means, variances and covariance calculated from the DNS data over a spatially filtered domain (5 cm filter width) in the formulations, and compare them to the corresponding diabatic contributions calculated directly from the DNS over the filtered domain. The filter width (5 cm) is slightly larger than the Taylor microscale (< 𝜆= 3.62-3.99 cm) in these simulations. The filter width was chosen to avoid significant impact of nonlinearity associated with molecular diffusion on the supersaturation while keeping the filter scale within the inertial subrange with application for LES in mind. The boundary points (100 grid points along all walls) are excluded in this comparison for simplicity.

Near the boundaries, the higher-order contributions become important because of the non-linearity, affecting the comparison. Figure 11 shows the results from this comparison. In the figure, each dot is from a volume of 5 3 cm 3 . The plots show that the modeled terms generally match the fully resolved DNS outputs. There is some deviation at lower magnitude points that could be the result of neglecting terms higher than second order in the simplified model. Larger (10 cm) and a smaller

(2 cm) filter widths were also tested (not shown). Our simplified model also performs well for these smaller and larger filter widths. However, the data scatter around the 1:1 line is more for the lowest filter width tested (2 cm < 𝜆), likely due to higher-order effects associated with molecular diffusion and/or higher statistical fluctuations associated with particle variability that are neglected in our model.

The DNS results presented in the previous section reasonably support foundational assumptions for stochastic subgrid-scale modeling (Gaussian PDFs of Lagrangian differences, exponential Lagrangian autocorrelation, delta-correlated Lagrangian differences, and Kolmogorov-Obukhov-Corrsin scaling in the inertial range). However, further evaluation of these assumptions is needed for higher Reynolds number turbulence. Criticism of the Lagrangian stochastic model commonly interprets supersaturation fluctuations in this model as a randomly fluctuating quantity without any structure (e.g., Prabhakaran et al. 2022). However, that is not entirely true. Lagrangian stochastic models (based on the Langevin approach) assume exponential Lagrangian correlation, which is supported by the DNS results in the previous section. The drift term in the model ensures exponential correlation with a characteristic timescale of the system.

Moreover, the Wiener forcing term is based on Kolmogorov-Obukhov-Corrsin scaling for scalars in the inertial range that is also consistent with the Lagrangian supersaturation fluctuations in our DNS.

The assumption of delta-correlated Lagrangian supersaturation gradients (not supersaturation fluctuations) is also reasonable, considering the very narrow correlation functions produced by the current DNS (not shown). Another common criticism of Lagrangian stochastic models is the assumption of Gaussian PDFs of supersaturation fluctuations, although it assumes Gaussian PDFs of Lagrangian supersaturation differences (not supersaturation itself). Consistent with past studies (e.g., Chandrakar et al. 2020;Thomas et al. 2021), the current DNS showed that supersaturation in the idealized adiabatic case (where water vapor and temperature are highly correlated) has a highly negatively skewed PDF. However, PDFs are closer to Gaussian for more realistic non-adiabatic cases where the water vapor and temperature correlations are weaker. The PDFs of Lagrangian supersaturation difference are closer to a Gaussian shape near the mode, but with tails significantly deviating from Gaussian owing to intermittency. Therefore, a Gaussian Wiener forcing for supersaturation fluctuations is a reasonable assumption for condensation growth, but may be less justified for the activation/deactivation process, where larger fluctuations in the PDF tails could be important.

To discuss implications of this work on subgrid scale modeling further, we write Lagrangian stochastic differential equations for subgrid vertical velocity, potential temperature, water vapor, and supersaturation fluctuations along a particle trajectory. The adiabatic expansion/compression contribution associated with vertical velocity fluctuations (not present in the current isobaric DNS) is also included in the temperature and supersaturation fluctuations for atmospheric cloud conditions. These equations can be expressed as (Paoli and Shariff 2009;Sardina et al. 2015;Chandrakar et al. 2016Chandrakar et al. , 2018;;Abade et al. 2018):

𝑑𝑇 ′ 𝑖 (𝑡) = -𝑇 ′ 𝑖 (𝑡) 𝜏 𝑡 𝑇 -𝐿/𝐶 𝑝 𝑠 ′ 𝑖 (𝑡) 𝜏 * 𝑐 𝑑𝑡 -𝑔/𝐶 𝑝 𝑤 ′ 𝑖 (𝑡)𝑑𝑡 + 2𝜎 2 𝑇 𝑑𝑡 𝜏 𝑡 𝑇 1/2 𝜉 𝑇 𝑖 (𝑡), (9) 𝑑𝑞 ′ 𝑣 𝑖 (𝑡) = -𝑞 ′ 𝑣 𝑖 (𝑡) 𝜏 𝑡 𝑞 + 𝑠 ′ 𝑖 (𝑡) 𝜏 * 𝑐 𝑑𝑡 + (1 -𝐶 2 𝑇,𝑞 𝑣 ) 2𝜎 2 𝑞 𝑣 𝑑𝑡 𝜏 𝑡 𝑞 1/2 𝜂 𝑞 𝑖 (𝑡) + 𝐶 𝑇,𝑞 𝑣 2𝜎 2 𝑞 𝑣 𝑑𝑡 𝜏 𝑡 𝑞 1/2 𝜉 𝑇 𝑖 (𝑡), (10) 𝑑𝑠 ′ 𝑖 (𝑡) = -𝑠 ′ 𝑖 (𝑡) 𝜏 𝑡 𝑠 -𝑠 ′ 𝑖 (𝑡) 𝜏 𝑐 + 𝛼(𝑇)𝑤 ′ 𝑖 (𝑡) 𝑑𝑡 + (1 -𝐶 2 𝑤,𝑠 ) 2𝜎 2 𝑠 𝑜 𝑑𝑡 𝜏 𝑡 𝑠 1/2 𝜉 𝑠 𝑖 (𝑡) + 𝐶 𝑤,𝑠 2𝜎 2 𝑠 𝑜 𝑑𝑡 𝜏 𝑡 𝑠 1/2 𝜁 𝑤 𝑖 (𝑡), (11) where, 𝑤 ′ 𝑖 , 𝑞 ′ 𝑣 𝑖 , 𝑇 ′ 𝑖 , and 𝑠 ′ 𝑖 are the subgrid vertical velocity, water vapor mixing ratio, temperature, and supersaturation fluctuations experienced by the 𝑖 𝑡ℎ particle along its trajectory; 𝜏 𝑡 is the Lagrangian correlation timescale of turbulent fluctuations (additional subscript indicates the velocity or scalar fields), 𝜏 * 𝑐 is a relaxation timescale due to the condensation as defined in Sec. 2b; 𝜎 2 𝑤 , 𝜎 2 𝑇 , 𝜎 2 𝑞 𝑣 , and 𝜎 2

𝑠 𝑜 are the strength of the Wiener terms for vertical velocity, temperature, water vapor, and supersaturation fluctuations; 𝛼 = 𝛽𝑔/(𝐶 𝑝 T) -𝑔/(𝑅 𝑑 T) is a constant coefficient that depends on temperature (𝑇); 𝜁 𝑤 𝑖 (𝑡), 𝜉 𝑇 𝑖 (𝑡), 𝜂 𝑞 𝑖 (𝑡), and 𝜉 𝑠 𝑖 (𝑡) are delta correlated independent Gaussian random noise functions with variance equal to unity. The delta-correlated noise terms assume a significantly shorter correlation timescale of Lagrangian supersaturation, water vapor, temperature, and velocity increments than the integral timescale and a Gaussian distribution of the increments within the inertial range. Independent noise functions for the Wiener forcing are used for different quantities since their time evolution is not perfectly correlated. However, 𝑞 ′ 𝑣 and 𝑇 ′ forcings are expected to be partially correlated in turbulent convection. Thus, the Wiener forcing in the 𝑞 ′ 𝑣 evolution has two orthogonal forcing terms, one with the same Gaussian noise as the 𝑇 ′ forcing and the other with an independent Gaussian noise term. The magnitudes of these two terms are based on the normalized covariance (𝐶 𝑇,𝑞 𝑣 = 𝑞 ′ 𝑣 𝑇 ′ 𝜎 -1 𝑞𝑣 𝜎 -1 𝑇 ) estimated from the subgrid model presented in Appendix A (Eqs. A1, A2, A3). As shown in Sec. 4b and Chandrakar et al. (2022), the covariance of water vapor and temperature is critical for supersaturation variability. Thus, the correlated and uncorrelated Wiener forcing terms for water vapor and temperature fluctuations are a key to representing supersaturation variability accurately. Similarly, the evolution of 𝑤 ′ and 𝑠 ′ fields are partly correlated through temperature variability driven by vertical velocity fluctuations. Prabhakaran et al. (2022) also discussed the importance of droplet transport correlated with growth for appropriate representation of the droplet size distribution evolution. Therefore, the 𝑠 ′ evolution equation also has two orthogonal Wiener forcing terms with magnitudes based on the normalized covariance (𝐶 𝑤,𝑠 = 𝑤 ′ 𝑠 ′ 𝜎 -1 𝑤 𝜎 -1 𝑠 ) estimated below using Eq. 13. For simplicity, the Wiener forcing term in the 𝑇 ′ equation is not separated into correlated and uncorrelated parts with the 𝑤 ′ equation. Nevertheless, the drift term associated with vertical velocity fluctuations ensures the partial correlation of temperature fluctuations with vertical velocity fluctuations at large scales.

There are two potential approaches to model the Lagrangian subgrid variability of supersaturation and droplet evaporation/growth: (a) using equations 8, 9, and 10, or (b) using equations 8 and 11. In the first approach, the properties of the independent scalars in turbulent flow are modeled explicitly without a need to deal with the properties of supersaturation separately. Thus, there is no need to assume a Gaussian PDF of Lagrangian supersaturation difference, determine the supersaturation autocorrelation timescale, or assume that the supersaturation follows the Lagrangian properties of independent scalars. However, this approach requires tracking an additional variable and solving an extra equation compared to the second approach. Thus, the first approach is more complicated and computationally expensive than the second. Moreover, temperature and water vapor PDFs may not be Gaussian in the regions influenced by the entrainment-mixing process. On the contrary, the supersaturation PDFs are closer to Gaussian for non-adiabatic cases (see Figs. 4 9). Thus, using the 𝑠 ′ equation with the second approach may have additional advantages. Further testing of these approaches by comparing to DNS statistics after implementation in an LES framework is left for future investigation.

As shown in the previous section, individual magnitudes of water vapor and temperature fluctuations do not change significantly before and after cloud formation. Thus, for the first approach, the drift term from condensation (second term) can be neglected, and the Wiener forcing terms for temperature and water vapor can be directly obtained from the Eulerian subgrid scalar variance and covariance model discussed in Appendix A (Eqs. A1, A2, A3):

and 𝐶 𝑇,𝑞 𝑣 ≈ ( 𝑞 𝑣 𝜃 -q𝑣 θ)𝜎 -1 𝜃 𝜎 -1 𝑞 𝑣 . The magnitude of the vertical velocity fluctuations (𝑤 ′2 = 𝜎 2 𝑤 ) can be estimated from a subgrid prognostic subgrid-scale turbulent kinetic energy model in LES. Note that the above relations assume that the Lagrangian subgrid scalar variances and covariances experienced by a particle at a given moment are approximately equal to their Eulerian subgrid magnitudes at the particle grid location.

For the second approach that uses direct modeling of the subgrid supersaturation, the strength of the Wiener forcing terms can be determined by a few additional steps discussed below. Using Ito calculus, equations for 𝑠 ′2 (= 𝜎 2 𝑠 ) and 𝑠 ′ 𝑤 ′ in the presence of phase change can be obtained using Eqs. 8 and 11:

and their quasi-steady-state solutions are:

Here, overlines indicate the ensemble average of different trajectories. As shown in the above equations, the condensation process and turbulent mixing of scalars decorrelate vertical velocity and supersaturation fields, but vertical velocity-driven adiabatic expansion/compression correlates the two. The strength of the Wiener forcing for supersaturation fluctuations (𝜎 2 𝑠 0 ) in Eq. 11 can be calculated by substituting 𝑠 ′2 ≈ ( 𝑠𝑠s s) in Eqs. 14 and 15, where ( 𝑠𝑠s s) is the Eulerian magnitude of supersaturation fluctuations in a grid cell obtained from Eq. 2. In turn, the terms in Eq. 2 for calculating ( 𝑠𝑠s s) can be obtained from the subgrid model given by Eqs. A1, A2, and A3. Note that the supersaturation autocorrelation timescale differs from the velocity autocorrelation timescale. As shown in Fig. 2, the supersaturation correlation timescale varies with external relative water vapor and temperature fluxes. Further investigation of the scaling of supersaturation autocorrelation timescale is needed for its appropriate parameterization as a function of water vapor and temperature fluctuation statistics and is left for future studies.

The representation of supersaturation variability in models is critical for subgrid interactions of cloud particles with turbulence fields. Using DNS of turbulent Rayleigh-Bénard convection as in the Pi Convection Cloud Chamber, we investigated different Lagrangian and Eulerian scalar and supersaturation statistics before and after cloud formation that are critical for subgrid modeling.

Lagrangian statistics of passive particles and cloud droplets (𝑅 > 5 𝜇m) were compared to study the impact of droplet gravitational settling. The effects of relative scalar fluxes (e.g., during the entrainment-mixing process) were investigated by changing the sidewall boundary conditions. A simplified model of the contributions from phase change to subgrid water vapor and temperature variances and covariance prediction was introduced and evaluated using DNS output.

The Lagrangian supersaturation autocorrelation displays an exponential shape similar to independent scalars. But its integral time is different than water vapor and temperature, and this integral time varies with the sidewall relative fluxes. Cloud formation suppresses the supersaturation correlation timescale but has little impact on the water vapor and temperature correlation timescales. The gravitational settling of droplets also decreases the supersaturation correlation timescale. The Lagrangian second-order structure functions of water vapor, temperature, and supersaturation are similar to that expected for a passive scalar in isotropic turbulence, although the Rayleigh-Bénard convection is non-isotropic (except in the convection core). An inertial range scaling of 𝐶𝜀 𝑥 𝜏 (or equivalently 𝐶 (2𝜎 2 𝑠 /𝜏 𝑡 𝑠 )𝜏 for supersaturation is evident. This is an important result given that such a scaling is a basis for Lagrangian stochastic subgrid models. Settling of droplets suppresses the magnitude of supersaturation fluctuations along droplet trajectories, but the shape of the Lagrangian structure functions remains nearly the same. PDFs of the Lagrangian supersaturation difference are also close to Gaussian near the mode, but the tails significantly deviate from Gaussian. This PDF shape is similar to past observations of the Lagrangian acceleration in isotropic turbulence (e.g., La Porta et al. 2001;Brown et al. 2009). These results therefore support using a Lagrangian stochastic formulation following the Langevin drift-diffusion approach for supersaturation fluctuations at the inertial range similar to that for velocity fluctuations but with vertical velocity-correlated and uncorrelated Wiener forcing terms.

Another quantity of interest in turbulent convection relevant for subgrid modeling is the shape of Lagrangian and Eulerian scalar PDFs. The shapes of the Lagrangian and Eulerian PDFs are not significantly different. The Lagrangian PDFs have slightly greater fluctuation magnitudes and extreme values than the Eulerian PDFs since particle trajectories include sampling of the regions near boundaries. The shape of the scalar and supersaturation PDFs is influenced by the sidewall relative scalar fluxes and cloud formation. Independent scalar PDFs (temperature and water vapor) in the convection core have nearly zero skewness in the Rayleigh-Bénard convection without sidewall fluxes (adiabatic case). Although close to Gaussian shape, the PDF tails deviate significantly from Gaussian (high kurtosis values). However, the PDF of supersaturation is negatively skewed and has a higher kurtosis value due to a nonlinear dependence on the individual scalars. The presence of negative external fluxes makes the water vapor PDFs positively skewed. Interestingly, the supersaturation PDFs have nearly zero skewness for these more realistic conditions in which external fluxes are present, contrary to the individual scalars. Moreover, the supersaturation PDFs are closer to Gaussian shape than the water vapor and temperature PDFs. For the non-adiabatic cases, the tails of the supersaturation PDFs deviate less from Gaussian, and the kurtosis is much less than that for water vapor and temperature fluctuations. Overall, these results suggest that assuming a Gaussian distribution of supersaturation fluctuations in subgrid condensation models is reasonably justified for non-adiabatic cases.

We also derived an expression for supersaturation skewness in terms of temperature and water vapor skewness and their third-order covariance. This analysis showed that the temperature and water vapor skewness have opposite contributions to supersaturation skewness. The third-order covariance terms are also significant and counteract one another. Cloud formation also changes the sign of the skewness from negative to positive for the adiabatic case and nearly zero to a positive value for saturated sidewalls. For the case with subsaturated sidewall boundaries, the PDF is closest to Gaussian shape.

The diabatic contributions through latent heating/cooling and water vapor changes in clouds from condensation/evaporation affect the magnitude of the scalar and supersaturation fluctuations. For water vapor and temperature variances and their covariance, the diabatic contributions depend on ( 𝑠𝑞 𝑣s q𝑣 ) and ( 𝑠𝜃s θ), which in turn depend on ( 𝑞 𝑣 𝑞 𝑣 -q𝑣 q𝑣 ), ( 𝜃𝜃θ θ), and ( 𝜃𝑞 𝑣θ q𝑣 ). The DNS results showed that cloud formation does not significantly impact water vapor and temperature fluctuations in the core of convection. However, the normalized covariance of water vapor and temperature is significantly impacted -increasing the covariance for non-adiabatic sidewalls and decreasing it for the adiabatic case. Thus, the magnitude of supersaturation variability decreases for the non-adiabatic cases and increases for the adiabatic case. The diabatic terms from condensation introduce a positive heat flux and negative water vapor flux at the phase relaxation timescale.

These fluxes decrease the covariance of water vapor and temperature fluctuations when they are highly correlated (and otherwise increase them), but increase the fluctuation magnitude of both scalars. A model accounting for these contributions in the scalar variance and covariance prognosis performed reasonably well when directly compared to DNS output.

This article provides detailed information on Lagrangian and Eulerian statistics of supersaturation critical for subgrid modeling of cloud microphysics. A major conclusion is that our DNS results support several key assumptions underpinning Lagrangian stochastic models of subgrid supersaturation variability, both in cloud and non-cloud conditions. This justifies the use of these models for representing supersaturation fluctuations and subgrid condensation in microphysics schemes, particularly Lagrangian particle-based schemes. The supersaturation fluctuations at subgrid scales of atmospheric models are significant for droplet activation and droplet size distribution broadening.

The proposed subgrid model will improve the representation of these turbulence effects in cloudresolving and LES models, especially during the turbulent entrainment-mixing process in clouds. Future work will focus on coupling the subgrid model based on the current study and Chandrakar et al. (2022) to Lagrangian microphysics in LES and evaluating it against droplet growth data from DNS of the Pi Cloud Chamber. Subsequently, we will apply the new subgrid modeling framework in simulations of atmospheric clouds. It is also critical to evaluate Lagrangian and Eulerian scalar and supersaturation statistics at higher Reynolds number turbulence to be further confident in the applicability of their scaling on a broader range of scales in atmospheric cloud conditions.

𝐷 ( 𝑞 𝑣 𝑞 𝑣 -q𝑣 2 ) 𝐷𝑡 = -1 ρ𝑎 𝜕 𝜕𝑥 𝑖 ρ𝑎 𝑢 𝑖 𝑞 𝑣 𝑞 𝑣 -ũ𝑖 𝑞 𝑣 𝑞 𝑣 + 2 q𝑣 ρ𝑎 𝜕 𝜕𝑥 𝑖 [ ρ𝑎 ( 𝑢 𝑖 𝑞 𝑣 -ũ𝑖 q𝑣 )] + 𝜕 𝜕𝑥 𝑖 𝐷 𝑣 𝜕 ( 𝑞 𝑣 𝑞 𝑣 -q𝑣 2 ) 𝜕𝑥 𝑖 -2𝐷 𝑣 𝜕𝑞 𝑣 𝜕𝑥 𝑖 𝜕𝑞 𝑣 𝜕𝑥 𝑖 -2𝐷 𝑣 𝜕 q𝑣 𝜕𝑥 𝑖 𝜕 q𝑣 𝜕𝑥 𝑖 -( 𝑞 𝑙 𝑞 𝑣 -˜ 𝑞 𝑙 q𝑣 ), (A2) 𝐷 ( 𝑞 𝑣 𝜃 -q𝑣 θ) 𝐷𝑡 = -𝜕 q𝑣 𝜕𝑥 𝑖 ( 𝑢 𝑖 𝜃 -ũ𝑖 θ) -𝜕 θ 𝜕𝑥 𝑖 ( 𝑢 𝑖 𝑞 𝑣 -ũ𝑖 q𝑣 ) -1 ρ𝑎 𝜕 𝜕𝑥 𝑖 ρ𝑎 ( 𝑢 𝑖 𝜃𝑞 𝑣 -ũ𝑖 θ q𝑣 -ũ𝑖 𝜃𝑞 𝑣θ 𝑢 𝑖 𝑞 𝑣 -q𝑣 𝑢 𝑖 𝜃) -2(𝐷 𝑣 + 𝐷 𝑡 ) 𝜕𝑞 𝑣 𝜕𝑥 𝑖 𝜕𝜃 𝜕𝑥 𝑖 -𝜕 q𝑣 𝜕𝑥 𝑖 𝜕 θ 𝜕𝑥 𝑖 + 𝐿 𝐶 𝑝 𝑞 𝑙 𝑞 𝑣 /Π -˜ 𝑞 𝑙 q𝑣 / Π -˜ 𝑞 𝑙 𝑞 𝑣 /Π -𝑞 𝑙 𝑞 𝑣 / Π -𝑞 𝑙 /Π q𝑣 -( 𝑞 𝑙 𝜃 -˜ 𝑞 𝑙 θ), (A3)

where 𝐷/𝐷𝑡 represent the material derivative, the dot represents the rate of change with time, 𝑢 𝑖 is a flow velocity component, 𝜌 𝑎 is the air density, 𝐷 𝑡 and 𝐷 𝑣 are molecular diffusivities of heat and moisture, 𝐶 𝑝 is the specific heat capacity at constant pressure, and 𝑞 𝑙 is the liquid water mixing ratio. The above budget equations are based on past studies (e.g., Wyngaard et al. 1978;Deardorff 1974a). The first term on the right-hand side is turbulent transport. The second is the turbulent production, the third is the molecular diffusion of variance, the fourth is the dissipation, and the last is the diabatic production from condensation and evaporation.