Sensitivity of Simulated Deep Convection to a Stochastic Ice Microphysics Framework

Ice microphysics parameterizations in models must make major simplifications relative to observations, typically employing empirical relationships to represent average functional properties of particles. However, previous studies have established that ice particle properties vary even in similar cloud types and thermodynamic environments, and it remains unclear how this so‐called “natural variability” impacts simulated deep convection. This uncertainty is addressed by implementing a stochastic framework into the Predicted Particle Properties microphysics scheme in the Weather Research and Forecasting model. The approach stochastically varies the coefficients of the mass‐size (m‐D) relationship (m=aDb) for unrimed and partially rimed ice. Using guidance from aircraft in situ measurements obtained during the Midlatitude Continental Convective Clouds Experiment (MC3E), the scheme samples from distributions of the prefactor (a) and the exponent (b) of the m‐D relationship. Simulations of two MC3E deep convective cases indicate that the stochastic m‐D scheme produces considerable variability of anvil cirrus cloud optical depth (τ) distributions, even for the same ice water path (IWP). Thus, the stochastic scheme produces variable cloud radiative forcing that is independent of IWP. This τ‐IWP relationship variability is nonexistent using the deterministic m‐D ensemble. Additional sensitivity tests are performed in which the fallspeed‐size relationship (V=cDd) is stochastically varied, resulting in variable precipitation amounts and rain rate distributions. Results are presented in the context of satellite and precipitation observations and include comparison with other ensemble configurations using perturbed initial and lateral boundary conditions and small‐amplitude noise added to the potential temperature field.


Introduction
Numerous studies have analyzed the sensitivity of simulated mesoscale convective systems (MCSs) to microphysics representation because these systems have a large impact on weather and climate (e.g., Bryan & Morrison, 2012;Chen & Cotton, 1988;Gilmore et al., 2004;Van Weverberg et al., 2013). Convective ice detrainment into the upper troposphere produces anvil cirrus shields that can significantly impact the atmosphere's radiation budget (Del Genio & Kovari, 2002;Fu et al., 1995;Machado & Rossow, 1993;Stephens, 2005;Zender & Kiehl, 1997). For example, Feng et al. (2011) showed that anvil cirrus shields from continental deep convective systems contributed up to 31% of the net longwave and shortwave (SW) cloud radiative forcing over the U.S. Southern Great Plains during several deep convection events. In addition to radiative impacts, Nesbitt et al. (2006) found that over 50% of precipitation in the tropics is attributable to MCSs and increases to 90% over certain geographical regions. The significant impacts of MCSs on the radiation budget, global annual rainfall, and the vertical transport of heat, momentum, moisture, and aerosols motivate a detailed understanding of MCS sensitivity to microphysical parameter variability.
Representing this variability is challenging because ice particle properties in cloud microphysics schemes are extremely simplified compared to the complex diversity of observed ice crystal habits and size distributions. For example, current bulk microphysical parameterizations typically relate ice particle mass, terminal fallspeed, and projected area to the maximum particle dimension through simple power laws with constant coefficients (e.g., Hong et al., 2004;Lin et al., 1983;Morrison et al., 2005;Milbrandt & Yau, 2005;Thompson et al., 2008). These schemes often use empirically derived values that are held fixed in time and space. However, many studies have shown that these parameters spatially and temporally vary even for similar thermodynamic and cloud conditions (e.g., Brown & Swann, 1997;Field et al., 2005;Heymsfield et al., 2004;McFarquhar & Black, 2004;McFarquhar et al., 2007). These findings suggest limitations of using constant parameters and that including "natural variability" of these parameters may be important for properly representing the nonlinear physics of cloud systems.
Previous microphysics sensitivity studies of MCSs have either compared simulations employing entirely different microphysics schemes or used a single microphysics scheme while changing individual parameter values that are constant in time and space. These approaches have been successful in elucidating the impacts of certain microphysical parameters on the structure and evolution of specific types of cloud systems. In particular, numerous studies have established the sensitivity of simulated convective systems to ice microphysical parameter settings (e.g., Adams-Selin et al., 2013;Bryan & Morrison, 2012;Gilmore et al., 2004;McCumber et al., 1991;Morrison & Milbrandt, 2011;Varble et al., 2011, among many others). However, these studies did not address the effects of ice parameter temporal or spatial variability within a simulation.
An intuitive method to account for the natural variability of parameter values is through a stochastic framework in which parameters vary in time and space. Previous stochastic modeling approaches have been well established in improving the representation of model uncertainty and ensemble spread relative to forecast error (e.g., Berner et al., 2017Berner et al., , 2009Jankov et al., 2017). Such approaches include stochastic perturbation of the model state variable physical tendencies (SPPT; Berner et al., 2015;Buizza et al., 1999;Palmer et al., 2009) and the stochastic kinetic energy backscatter scheme (Berner et al., 2009(Berner et al., , 2011Mason & Thomson, 1992;Shutts, 2005), which uses a stochastic framework to account for upscale-propagating errors caused by unresolved subgrid-scale processes.
A criticism of these approaches, particularly SPPT, is the lack of a direct connection to uncertainty in the parameterized subgrid-scale physics. Thus, stochastically perturbing parameters directly within parameterization schemes has also been explored (e.g., Berner et al., 2015;Christensen et al., 2015;Jankov et al., 2017). However, only a few studies have analyzed the impacts of including stochastic microphysics. McFarquhar et al. (2003) applied a stochastic scheme in a single-column model for sampling ice effective radius and found an upscale impact on cloud radiative forcing. In Kuell and Bott (2014), the initial ice nuclei concentration and the turbulent collection kernel for rain formation were varied stochastically, which contributed (among other stochastic approaches) to increasing precipitation forecast skill for two simulations. Qiao et al. (2018) analyzed idealized supercell simulations with a scheme that stochastically perturbed the microphysics heating tendency and another that perturbed the particle size distribution (PSD) intercept parameter of a single-moment bulk scheme and found improved forecasts of updraft helicity intensity. Our study addresses impacts on simulated MCSs by stochastically varying microphysical parameters that are typically held constant in bulk schemes, but in contrast to previous studies the parameter perturbations are guided by observations. The purpose of this study is to characterize the sensitivity of cloud and precipitation structure to microphysical parameter perturbations in deep convection simulations. To this end, we address two separate but related aspects: (1) the effects of systematic model error, where a parameter is constant but set to a biased value; and (2) random model error, where a parameter exhibits spatiotemporal variations due to spatiotemporal inhomogeneities and uncertain flow dependence. A unique aspect of this study is to test the sensitivity to both random and systematic parameter uncertainties. Unlike most previous studies, we also robustly characterize the significance of sensitivities by comparison to an ensemble employing only small-scale, small-amplitude perturbations to the potential temperature field (i.e., noise). Ensemble simulations are then formed with parameter values stochastically varied in time and space to address "natural" parameter variability or with parameter perturbations that are held fixed in time and space within a simulation to address the role of parameter bias. For microphysical parameter variability context, we also include an ensemble with perturbations to large-scale initial and lateral boundary conditions. An ensemble framework is necessary here to explore the range of realizations associated with parameter perturbations. We emphasize that this application is different from typical ensemble evaluation focused on forecast spread relative to error.
We also compare the simulations to observations to determine if microphysical parameter perturbations lead to more consistent comparisons with the mean behavior of the ensemble relative to the baseline fixedparameter simulation. Characterizing the spread from different stochastic realizations can also be informative as to how these microphysical variations impact broader system behavior and whether this spread encompasses observed values. As such, our analysis is primarily restricted to fields that can be directly compared to observations (viz., precipitation and satellite-derived cloud products).
While there are dozens of parameters in a given microphysics scheme that could be stochastically perturbed, we investigate only a few. Recent supporting work by Finlon et al. (2018) has characterized the variability of the mass-size (m-D) relation using in situ observational retrievals. Our current study therefore focuses primarily on varying the coefficients of the m-D relationship wherein the sampled parameter values are guided by observed distributions. Several studies have also suggested the importance of the ice fallspeed-size relation (V-D) in simulating convective systems (e.g., Bryan & Morrison, 2012;Lord et al., 1984;McCumber et al., 1991;Morrison & Milbrandt, 2011). For example, Forbes and Clark (2003) tested the sensitivity to halving and doubling the ice fallspeed in a mesoscale model and found that frontal rain band development was modified by associated alterations in sublimational cooling. The sensitivity of simulations to perturbing this relation is also investigated.
The remainder of the paper is organized as follows. The scheme development methodology, including the observational basis for stochastically perturbing the m-D relationship parameters, is described in section 2. Observational and modeling data sets are discussed in section 3. The experimental design is presented in section 4. Finally, results are presented in section 5 with a discussion and conclusions in section 6.

P3 Microphysics Scheme
The stochastic framework is implemented within the Predicted Particle Properties (P3) bulk microphysics scheme (Milbrandt & Morrison, 2016;Morrison & Milbrandt, 2015). The P3 scheme includes one or more "free" ice categories; here we use a single ice category. For each category, four prognostic ice variables, the total ice mass, rime mass, rime volume, and number mixing ratios, allow evolution of various properties such as bulk density and fallspeed. The mass-size (m-D) relationship varies smoothly for up to four different size ranges of the PSD: (1) small, spherical ice, (2) nonspherical, unrimed ice (grown by vapor diffusion and/or aggregation), (3) partially rimed ice, and (4) fully rimed ice (i.e., graupel or hail). We refer to Figure 1 of Morrison and Milbrandt (2015) or Figure 1 of Dietlicher et al. (2018) for graphical context of how the m-D relationship varies for each PSD size range. For each range, all m-D relationships follow a power law: For spherical particles, the parameter a is proportional to the bulk density of ice (ρ), defined as particle mass divided by the volume of a sphere for a given D, and is given by Small ice crystals are assumed to be spherical (i.e., b=3) with an effective density equal to the bulk density of solid ice (917 kg/m 3 ). Graupel/hail is also assumed to be spherical with b=3, but the effective density varies and the a parameter is completely determined by the prognostic ice variables. On the other hand, a and b are determined empirically for nonspherical, unrimed ice, and the particle dimension is defined as the maximum particle dimension. In the standard version of P3, these a and b parameters are specified as from the relationship in Brown and Francis (1995) (b = 1.9, a = 0.0185 kg/m b ), modified by Hogan et al. (2012) (b = 1.9, a = 0.0121 kg/m b ). The deterministic scheme here uses slightly different parameter values (b = 2.1 and a = 0.0386 kg/m b ), to be consistent with the mean values of a and b in the stochastic scheme. For partially rimed ice, with an effective density between that of nonspherical, unrimed ice and fully rimed (graupel/hail) ice particles for a given D, b is equal to that for nonspherical, unrimed particles, while a is increased by a factor that depends on the rime mass fraction (equal to the ratio of the rime mass and total ice mass mixing ratios). Physically, this corresponds to assumptions that rime mass accumulates along the minor axis of ice particles, and the rime mass fraction is constant with D (see Morrison and Milbrandt 2015 for details). If the riming is heavy and a increases to the point that the mass of a partially rimed ice particle is equal to that of graupel/hail for a given D, then the ice is considered "fully rimed" and further growth follows the m-D relationship for graupel/hail.
Here we stochastically perturb the empirical a and b parameters, and thus, we modify the m-Drelationships for nonspherical, unrimed ice and partially rimed ice. The m-D relationships for small spherical ice and graupel/hail are not directly modified. We note that perturbing the a and b parameters in P3 requires modifications to associated lookup table creation code, and thus, the version of the scheme used herein is different from that in the default Weather Research and Forecasting (WRF) code.
Ice fallspeed (V) in P3 follows from Mitchell and Heymsfield (2005), such that V-D is related to the ratio of particle mass to projected area. The particle mass follows from the m-D relationship described above, whereas the projected area (A) follows from a power law A-D relationship similar to the m-D relationship (see Morrison and Milbrandt, 2015). Thus, perturbing the m-D relationship automatically perturbs the ice fallspeed of unrimed and partially rimed ice. However, we test additional sensitivity to ice fallspeed by applying a stochastic multiplicative factor to the V-D relationship for all ice. As such, the stochastic fallspeed scheme applies to rimed ice as well. Wu and McFarquhar (2016) describe the instrumentation used to collect the cloud microphysical measurements and the manner in which the data were processed. Note that the diameter definition used for in situ data is the maximum particle dimension, D max , consistent with the P3 scheme. Following that study, the two-dimensional cloud probe (2D-C; 125 μm <D max < 1 mm) and high-volume precipitation sampler (HVPS; D max > 1 mm) were combined to give a composite size distribution, with the probes agreeing within an average of 5% in the overlap range. The 2D-C had antishattering tips (Korolev et al., 2011) and algorithms based on the interarrival times of particles in the sample volume were applied to remove remaining artifacts following Field et al. (2003Field et al. ( , 2006. Data were processed using the University of Illinois/Oklahoma Optical Array Probe Processing Software (McFarquhar et al., 2017 with particle habits identified using a modified Holroyd III (1987) identification scheme; see Wu and McFarquhar (2016) for details. A deep-cone Nevzorov total water content sensor (Korolev et al., 1998) provided the total water content, which is equivalent to the ice water content for ice clouds. Although previous studies have shown that this probe can underestimate the ice water content, especially in the presence of larger particles due to them bouncing out of the cone (e.g., Wang et al., 2015), it does provide some information to constrain bulk mass and reduces the uncertainties in estimating a three-dimensional volume from a two-dimensional particle.

In Situ Observational Guidance
Measurements were obtained in the stratiform region of an MCS that occurred on 20 May 2011 during the Midlatitude Continental Convective Clouds Experiment (MC3E) field campaign (see section 3.1 for details on the campaign). Note that the (stratiform) sampling region is one reason for limiting perturbations to only unrimed and partially rimed ice. By doing this, we limit representative issues that would arise by applying perturbations to ice formed in the convective regions. Moreover, the a-b coefficients for rimed ice already vary in P3 based on the predicted rime density. That said, some representative issues may still exist due to sampling limitations, such as using observations from only one case (20 May), and thus, future work should attempt to aggregate samples from many cases.
The data processed at 1-s resolution are used in the analysis; although 1-s in situ PSDs are not ideal because of statistical sampling issues, they are used in order to increase information available about the small-scale variability within the cloud. To derive a-b coefficients for each 1 s of data, a b value within the range of 1.5 to 2.7 is first specified, with this range encompassing the range of most previously published b values . The unique a corresponding to each b is then determined by forcing the mass integrated from the PSD and the m=aD b relation to match that measured by the Nevzorov probe. The result is a distribution of a-b pairs that served as the sampling distribution for the stochastic scheme. For each a and b pair, an effective density can be derived as a function of maximum dimension using the volume of an equivalent spherical particle (with diameter D eq ) from the equation: Note that because a and b are strongly correlated (e.g., Finlon et al., 2018), their values cannot be sampled independently and the appropriate pairing must be used to determine the effective density. Thus, we instead sample from a distribution of effective densities at a specific size (here, 500 μm is chosen). Figure Note that equation (5) using the ρ 500μm value at b = 2.1 is used for the nonstochastic (deterministic a-b) simulations as opposed to the standard a parameter used in the default P3 scheme.

Stochastic Perturbation Scheme
A detailed description of the scheme used for stochastic sampling is provided by Berner et al. (2015) and Jankov et al. (2017), so only a brief description is provided here. For a stochastically varying parameter γ, values are obtained following where γ is the mean and r(x,y,t) is the perturbation field in spectral space defined by time (t) and a twodimensional Cartesian plane (x,y). The two-dimensional perturbation field r(x,y,t) is given by in which k and l are the (K + 1)(L + 1) wavenumber components in the x and y directions, respectively, r k,l is the spectral coefficient, and X and Y represent the number of grid points in the x and y directions, respectively. The spectral coefficient r k,l evolves over time as a first-order autoregressive process. Explicit description of this equation is given in Equation A4 of Berner et al. (2015). For brevity, we note that r k,l contains three terms: (1) a linear autoregressive parameter that is a function of the model time step and a prescribed temporal decorrelation time, (2) a wavenumber dependent noise amplitude that is a function of a prescribed length scale and standard deviation, and (3) a Gaussian white-noise process. Therefore, prescribing a temporal decorrelation time and a length scale allows the stochastic perturbations to evolve in time and space with an assigned degree of "memory." The result is a stochastically sampled, two-dimensional, time-varying field of correlated parameter values that follow a Gaussian distribution with a prescribed mean and standard deviation. Following this approach, fields of stochastically varying b and ρ 500μm for the m-D relationship are generated, using the means and standard deviations given in section 2.2. Subsequently, a is calculated from the b and ρ 500μm values following equation (5). For simplicity, parameter values are only varied horizontally; vertical variability is not considered.
For the stochastic ice fallspeed scheme, ice fallspeed is varied by applying a stochastic multiplicative factor F to the fallspeed-size (V-D) relationship. Values are sampled from a Gaussian distribution in log 2 space to ensure consistency with the mean V-D relationship. Thus, the mean of log 2 (F) is 0 and the standard deviation is specified to be 0.5, meaning that the range of F with ±2 standard deviations from the mean lies between 0.5 and 2. Although the approach for representing stochastic variability in the ice fallspeed relationship is more ad hoc compared to the ice m-D relationship, observational studies presented in the literature are used as rough guidance. Laboratory studies presented in Locatelli and Hobbs (1974) showed that the fallspeeds of ice crystals and their assemblages can vary by a factor of about 2 for a given D. Using in situ aircraft observations and a theoretical calculation of ice fallspeed, Heymsfield et al. (2007) showed that the fallspeed of cirrus ice crystals range from about 0.5 to 2.5 m/s. Multiple relationships from the literature were plotted by Rosenow et al. (2014) to show this variation. In the context of these studies, a multiplicative factor with ±2 standard deviations from the mean, between 0.5 and 2, is within reasonable physical bounds.

MC3E
Simulations are compared with observations from MC3E (Jensen et al., 2016) collected during the spring of 2011. MC3E was a joint field campaign held by the U.S. Department of Energy's Atmospheric Radiation Measurement Program and the National Aeronautics and Space Administration's Global Precipitation Measurement mission in north central Oklahoma. This study focuses on two MC3E case studies of contrasting synoptic-scale forcing and convective mode.
The first case was a squall line that occurred on 20 May, in which the forcing mechanism was a cold front associated with a strong upper-level trough usurping a dry line in western Oklahoma. This system organized around 0600 UTC and moved eastward across the state for the next 15 hr, producing a distinct leading convective line and area of high reflectivity, a wide trailing stratiform region, and an expansive anvil cirrus shield covering most of Oklahoma. The simulations were initialized at 0000 UTC and run through 2100 UTC on 20 May.
The second case occurred from 23-25 May and was fundamentally different than the 20 May case. Soundings taken at six locations across Oklahoma show significant low-level helicity, 25 m/s 0-to 6-km vertical wind shear, and convective available potential energy values exceeding 4,000 J/kg for some times and locations.
The key feature for deep convection initiation was the presence of a pronounced dryline over western Oklahoma that moved eastward. Supercells initiated along the dryline beginning around 2000 UTC on 23 May. By 0000 UTC on 24 May, the cells in northern Oklahoma began to organize into a well-defined bow echo that drove southeastward over the next 6 hr. Convection then largely ceased for the next 12 hr. At around 1900 UTC on 24 May, convection initiated in western Oklahoma along a north-south line associated with a strong closed upper-level low and associated midlatitude cyclone. By 2100 UTC, a squall line with a forward anvil had formed and moved eastward across the state over the next 9 hr. Simulations of this event are run from 0000 UTC on 23 May to 0900 UTC on 25 May.

Observational Data Sets
Hourly rainfall estimates are provided by the Arkansas-Red Basin River Forecast Center (ABRFC; ARM, 1994) on a 4-km grid. The entire ABRFC domain is shown by the orange polygon in Figure 2. ABRFC precipitation estimates combine WSR-88D NEXRAD radar precipitation estimates with rain gauge reports and have been used in a number of precipitation analyses over the Southern Great Plains (Sun & Krueger, 2012;Tang et al., 2017;Varble, 2018;Xie et al., 2014;Zhang & Klein, 2013). Additional rainfall observations are from the gauge-corrected NEXRAD National Mosaic and Multisensor Quantitative Precipitation 10.1029/2019MS001730

Journal of Advances in Modeling Earth Systems
Estimate (hereafter Q2; Zhang et al., 2011), available at 1-km resolution but degraded here to 4-km resolution for comparison with ABRFC. The Q2 domain spans much of the conterminous United States, but only the portion within the ABRFC domain is considered here. Sensitivity tests were performed with precipitation analysis on the larger Q2 domain, but results were qualitatively similar to the ABRFC-domain-limited analysis in the context of model comparison.
The Satellite Cloud and Radiation Property retrieval System (SatCORPS) is used for satellite comparisons (ARM, 2011). The SatCORPS data sets use multiple algorithms to provide cloud property information from meteorological satellite data (Minnis et al., 2008). The domain used for SatCORPs measurements is shown by the blue box in Figure 2. For anvil detection (see section 5.2), SatCORPS data are used in combination with NEXRAD reflectivity, which is gridded using the Python Arm Radar Toolkit (Helmus & Collis, 2016) over the SatCORPS domain to comparable resolution.
Finally, broadband SW radiometer observations of downwelling SW irradiance and estimated clear-sky downwelling SW irradiance at the Southern Great Plains supplemental sites (red triangles in Figure 2) are used to compute cloud radiative forcing (ARM, 1997). Clear-sky downwelling SW irradiance is derived by applying a fitting algorithm to the SW radiometers (Long & Ackerman, 2000).

Model Setup
The Advanced Research Weather Research and Forecasting model Version 3.8 (Skamarock et al., 2008) is used for all simulations. WRF-Advanced Research Weather is a fully compressible, nonhydrostatic mesoscale model and is run at 3-km horizontal grid spacing with 51 vertical levels for the 3,000-km × 2,700-km domain shown in Figure 2 (black polygon). Vertical levels use an eta coordinate system in which layer depths increase with height in order to better resolve boundary layer processes. Aside from the ensemble employing large-scale initial and boundary condition perturbations (see section 4 for details), all simulations are initialized with a control member (CNTL) of the Second-Generation Global Medium-Range Ensemble Reforecast Dataset (GEFSR; Hamill et al., 2013). The GEFSR data set uses approximately 0.5°grid spacing and 42 vertical levels for 11 members (10 perturbed members + CNTL) and is run each day from 0000 UTC initial conditions. The physics parameterizations common to all simulations include the Mellor-Yamada-Janjic boundary layer scheme (Janjić, 1994), the Unified Noah Land Surface Model (Chen & Dudhia, 2001), and the Rapid Radiative Transfer Model for General circulation models longwave and SW radiation schemes (Iacono et al., 2008). The Mellor-Yamada-Janjic scheme performs vertical diffusion. Horizontal diffusion is calculated via a 2-D Smagorinsky scheme in which the diffusion coefficient is proportional to the product of deformation magnitude and the square of the turbulent mixing length (see Skamarock et al., 2008;Xue et al., 2000).

Experimental Design
A summary of the simulation ensembles is given in Table 1. Fifteen simulations are performed for each of the two stochastic microphysics configurations described in section 2 (stochastic a-b, hereafter STOCH-AB, and stochastic ice fallspeed, hereafter STOCH-VI). Inclusive of these 15 members for each scheme are 5 members per spatial (L) and temporal (t) autocorrelation scale: (1) t = 600 s, L = 18 km (appendage -SHORT), (2) t = 3,600 s, L = 100 km (appendage -MID), and (3) t = 10,000 s, L = 300 km (appendage -LONG).
Because there is currently no observational evidence for spatiotemporal autocorrelation scales of the m-D parameters, we test the scale sensitivity by employing this wide range of scales. We note also that using five ensemble members per stochastic ensemble may result in undersampling of some metrics (in particular, spread as quantified by ensemble standard deviation). However, a test was performed in which 15 additional members were run for the STOCH-AB-LONG ensemble (20 members total). For the metrics presented here, these results showed that only five to eight members were needed for the spread to saturate with respect to ensemble size and that using more members would not qualitatively change our results.
Stochastic simulations are compared to an ensemble with five members employing only small-amplitude grid-scale perturbations to the initial potential temperature (θ) field in the bottom seven model layers (roughly representative of the boundary layer). Perturbations to θ are applied using a random-number generator to create a Gaussian distribution centered at 0 with 2σ ≈ 0.1 K, where σ is the standard deviation. This ensemble (hereafter θ-pert) is a way to quantify the spread solely from growth (both in terms of scale and amplitude) of initially small-scale, small-amplitude errors (growth which has been shown to saturate at the convective scale after approximately 6 hr, e.g., Wang et al., 2012). We consider one measure of meaningful sensitivity from stochastic microphysics by whether or not there is a signal in the analyzed output that is larger than that produced by small-amplitude, grid-scale noise added to the initial θ field.
In addition to θ-pert, an ensemble with large-scale perturbations to the initial and lateral boundary conditions is created. This perturbed initial and lateral boundary condition ensemble (hereafter ICBC) is generated using the nonstochastic P3 scheme but initializing (and applying boundary conditions every 3 hr) with the 10 GEFSR members.
The fixed-parameter ensembles employ values at the lower and upper parameter ranges and thus are expected to produce large sensitivities. Specifically, the FIXED-AB b parameters employ ±2 standard deviations from the mean b, and the associated ρ 500μm values range from ±1-2.5 standard deviations from their respective means for a given b value. Although using these extreme ranges has been the typical approach for examining microphysical parameter sensitivity, it does not reflect natural parameter variability that is expected to be Gaussian like, for example, as seen for the observationally derived ρ 500μm parameter distribution discussed in section 2.2. To address this point, two final ensembles are created. In the first (FIXED-AB-

Journal of Advances in Modeling Earth Systems
SAMPLED), five b−ρ 500μm pairs are randomly sampled from the observationally constrained a and b distributions, which are thus representative of a random sample from the distribution. Similarly, the FIXED-VI-SAMPLED ensemble contains five members that are randomly sampled from the distribution of ice fallspeed multiplicative factors. These randomly sampled fixed-parameter ensembles are analogous to the stochastic ensembles in the limit of a very large spatiotemporal autocorrelation scale relative to the size of the domain and model integration time.
We note that the variability produced by the fixed-parameter ensembles is expected to be larger than the stochastic simulations since the stochastic scheme samples from a Gaussian distribution with perturbations from the mean, while each fixed-parameter ensemble member uses constant values not equal to the Gaussian mean. As such, the fixed-parameter ensemble members may be expected to yield systematic changes in the chosen metrics. In addition, variability of the ICBC ensemble is expected to be (much) larger than that produced by the stochastic scheme since ICBC perturbs the large-scale thermodynamic and kinematic environment, while the stochastic scheme only acts on microphysical processes where ice hydrometeors are present. A comparison of the stochastic ensembles with ICBC and fixed-parameter ensembles is performed only to add context from these other commonly used ensemble configurations. We also expect the sensitivity produced by the FIXED-AB-SAMPLED and FIXED-VI-SAMPLED ensembles to be smaller than that from the FIXED-AB and FIXED-VI, respectively, since the sampled parameter ensembles use values that are sampled from a Gaussian distribution whereas FIXED-AB and FIXED-VI use values at the bounds of physically reasonable parameter ranges.

Precipitation Structure 5.1.1. Control Simulations of 20 and 23-25 May 2011 Events
Simulations of the 20 and 23-25 May 2011 events produced overall system evolution similar to observed, but with precipitation field biases regardless of the model setup. Figure 3 shows the ABRFC observed (Figures 3a,  3c, 3e, and 3g) and simulated (Figures 3b, 3d, 3f, and 3h) hourly accumulated precipitation at specific times for each event. It also shows the observed (Figures 3i and 3k) and simulated (Figures 3j and 3l) temporal evolution of volumetric precipitation via Hovmöller diagrams, where volumetric precipitation is calculated as the total hourly accumulated precipitation for a grid cell multiplied by the area of the grid cell (units of millimeters per square kilometer per hour) and is summed over the meridional extent of the analysis subdomain for each longitude. Simulated results are shown for the nonstochastic control (hereafter, CNTL) member, which uses the nonstochastic P3 scheme and a single member of the GEFSR ensemble for initial and boundary conditions that was chosen as the control member because it gave precipitation results closest to observations.
At 1000 UTC on 20 May, there is less simulated hourly accumulated precipitation compared to what was observed, and the simulated precipitation field is displaced further north. By 1600 UTC on 20 May (Figures 3c and 3d), the simulation produces more precipitation ahead of the previous leading convective line than was observed, while the remnant leading convective line is more disorganized and produces weaker precipitation. Hovmöller diagrams (Figures 3i and 3k) show the low bias in precipitation along the dominant convective line before 1400 UTC, with a spatial separation of the precipitation field after this time into two distinct maxima. In contrast, the observed precipitation field is more spatially continuous throughout the time period.
The CNTL simulation also captures general features of the observed 23-25 May event. These include two precipitation periods before 1800 UTC on 23 May, followed by the initiation of supercells around 2100 UTC, and the development of a long-lived bow echo between 0000 and 1200 UTC on 24 May. Figures 3e and 3f show the observed and simulated 0300-0400 UTC accumulated precipitation associated with the bow echo. The simulated system produces less rainfall, a less expansive precipitation field, and a northwestward displacement compared to what was observed. By 2100 UTC on 24 May, a squall line developed and moved eastward across the analysis domain until 0900 UTC on 25 May. Hovmöller plots of the meridionally summed precipitation (Figures 3k and 3l) show that the model captures the general timing of this event but produces much less precipitation than observed, never fully achieving the observed linear squall line structure (Figures 3g and 3h). We note that tests were performed with initial and boundary conditions from the National Centers for Environmental Prediction Final Operational Global Analysis data set, but similar biases were produced.

10.1029/2019MS001730
Journal of Advances in Modeling Earth Systems

Fixed-Parameter Ensembles Compared to ICBC and θ-pert Ensembles
Domain-accumulated volumetric precipitation is shown as a function of time in Figure 4 for the ICBC, θ-pert, FIXED-AB, and FIXED-AB-SAMPLED ensembles. Both the 20 and 23-25 May cases are shown. All 20 May ICBC members underproduce precipitation, although they capture the temporal evolution of precipitation fairly well (Figure 4a). This is evident for the 23-25 May case as well (Figure 4c). The ICBC ensemble envelops the CNTL member for 20 May with a factor of 3 maximum spread at 1800 UTC. In contrast, the θ-pert ensemble has little deviation from the CNTL member for 20 May (Figure 4b).
The variability in accumulated precipitation for the 20 May FIXED-AB ensemble (Figure 4e) can partly be explained in terms of the different m-D relationships. For context, these m-D and ρ-D relationships are shown in Figures 5a and 5d. Figure 4e shows that FIXED-AB members with b = 2.7 produce more

Journal of Advances in Modeling Earth Systems
precipitation than members with b = 1.5. For a given b parameter, members with larger ρ 500μm produce larger precipitation amount as well. Simulations with b = 2.7 distribute more mass at larger particle sizes ( Figure 5a) compared to b = 1.5, and this larger m-D slope for b = 2.7 likely leads to greater size sorting, increased collision-coalescence of precipitation-sized particles, increased fallspeeds near the melting level, and thus greater downward condensate mass flux and precipitation accumulation.
The FIXED-AB-SAMPLED ensemble for 20 May is also shown (Figure 4f) with the respective m-D and ρ-D relationships in Figures 5b and 5e. This ensemble produces less spread than FIXED-AB, as expected given the smaller range of b and ρ 500μm parameters. Despite this smaller spread, the m-D and ρ-D relationships show that, similar to the FIXED-AB ensemble, the members with larger b and ρ 500μm lead to larger amounts of precipitation.
Precipitation evolution for the 20 May case was strongly constrained by significant synoptic forcing such that differences among the FIXED-AB ensemble members can at least be partially explained by precipitation microphysical processes directly responding to a-b perturbations. However, establishing the same trends in precipitation structure for the 23-25 May case is more challenging given the weaker synoptic forcing and the nonlinear, indirect impact of dynamical-microphysical feedbacks on surface precipitation. For example, note that the b = 1.5 simulations produce more precipitation between 0600 and 1200 UTC on 24 May (Figure 4g) compared to the b = 2.7 simulations (opposite of the 20 May FIXED-AB results). Spatial precipitation analysis (not shown) suggests that differences in precipitation accumulation at this time between the b = 1.5 and b = 2.7 simulations result from differences in bow echo evolution (the main source of precipitation at this time). However, these bow echo evolution differences cannot robustly be attributed to the varying b parameters since Figure 5d also shows one θ-pert ensemble member that behaves similarly to the b = 2.7 simulations.

Stochastic a-b Parameter Ensembles and Dependence on Spatiotemporal Autocorrelation Scale
The stochastic P3 scheme produces parameter values based on the Fourier modes in equation (7), as illustrated in Figure 6. Note that the ρ 500μm and b parameters are calculated everywhere based on the smooth Fourier modes but are only used in the scheme if ice is present at a given time and grid location. Figure 6 shows examples of the ρ 500μm and bparameter fields at 1200 UTC on 20 May for individual ensemble members with the short and long autocorrelation scales. The longer autocorrelation scale members also have parameter fields that evolve more slowly over time.
Results from the fixed-parameter ensembles show that simulated precipitation is slightly sensitive to the choice of a-b parameters. The stochastic ensembles produce precipitation accumulation variability as well, but less variability compared to the fixed-parameter ensembles, as expected (Figure 7). Precipitation accumulation from the stochastic ensembles for 20 May deviates only slightly from the CNTL member. A larger relative spread is seen at the peak of the precipitating period for 23-25 May, but recall that even θ-pert produces some spread at this time ( Figure 4d) and this is likely due to changes in the dynamical evolution of the bow echo.

Journal of Advances in Modeling Earth Systems
Spread comparisons, as quantified by standard deviation (σ), among all of the ensembles are quantified in Figure 8. The stochastic a-b ensembles produce slightly larger σ than θ-pert for most times (Figures 8b  and 8e), with STOCH-AB-MID and STOCH-AB-LONG producing slightly larger σ than STOCH-AB-SHORT. The normalized spread (ensemble σ divided by the mean) for the 20 May case produced by the stochastic schemes ranges from 2-8%, while θ-pert remains below 3% for most of the time period (Figure 8c). Despite larger σ for the stochastic ensembles compared to θ-pert, there is little variability in the ensemble means (Figure 8a), suggesting that the mean behavior for this case is more strongly controlled by the synoptic-scale forcing and that microphysics perturbations employed here do not strongly impact precipitation evolution.
Some similar trends in σ exist for the 23-25 May case compared to the 20 May case (Figures 8d-8f). For example, the STOCH-AB-MID and STOCH-AB-LONG ensembles generally produce larger σ compared to STOCH-AB-SHORT. However, the θ-pert ensemble for this case produces larger σ than STOCH-AB-SHORT at the peak of the precipitating period (0600 UTC on 24 May). Normalized σ ranges from 10-55% for the stochastic ensembles in this case, with larger normalized σ occurring during periods of lower precipitation ( Figure 8f). Fixed-parameter ensembles also clearly produce larger σ compared to the stochastic ensembles for both the 20 and 23-25 May cases, which is expected.
Some of the domain-accumulated precipitation variability shown in Figure 4 from the fixed-parameter ensembles appears to be correlated with variability in hourly rain rate distributions. Cumulative distribution functions (CDFs) of rain rate contributions to total rainfall (defined as the cumulative sum of hourly rain rates normalized by the total sum of precipitation) are shown in Figure 9 for three separate time periods: (1) the 20 May squall line from 0700-2100 UTC (Figures 9a-9d), (2) the supercells and bow echo that occurred from approximately 1800 UTC on 23 May until 0900 UTC on 24 May (Figures 9e-9h), and (3) the midlatitude cyclone and associated squall line that moved through the analysis domain from approximately 1800 UTC on 24 May until 0900 UTC on 25 May (Figures 9i-9l). The 20 May FIXED-AB

Journal of Advances in Modeling Earth Systems
ensemble member with b = 1.5 and ρ 500μm = 50 kg/m 3 (red line) produced the smallest hourly rain rate contributions to total rainfall ( Figure 9c) and also produced the smallest amount of domain-accumulated precipitation in Figure 4. However, there is little variability otherwise in the FIXED-AB and FIXED-AB-SAMPLED ensembles for the 20 May case. There is also very little variability in hourly rate contribution CDFs for the supercells/bow echo period (Figures 9e-9h) in the fixed-parameter ensembles, while variability is large for the ICBC ensemble ( Figure 9e). During the midlatitude cyclone/squall line time period, there is some correlation between greater high rain rate contributions and greater accumulated precipitation (yellow line, Figure 9k) and greater low rain rate contributions being correlated with lesser accumulated precipitation (red line, Figure 9k).
Rain rate distribution variability for the stochastic ensembles is sometimes very small and sometimes as large or larger than that for the fixed-parameter ensembles ( Figure 10). For the 20 May case (Figures 10a-10c), STOCH-AB-SHORT produces more variability in rain rates than STOCH-AB-MID or STOCH-AB-LONG, but the spread is very small for all autocorrelation scales. During the supercell/bow echo event (Figures 10d-10f), there is little variability in hourly rain rate contribution CDFs, likely because the precipitation during this case results primarily from convective processes (riming) as opposed to stratiform processes (vapor deposition). Similar variability is produced among all stochastic ensembles for the 24-25 May midlatitude cyclone/squall line event, with median rain rates varying by 2 mm/hr, similar to FIXED-AB for this event (Figure 9k). However, median rain rate variability for the stochastic ensembles in this case is comparable to that produced by θ-pert (cf. Figure 9j and Figures 10g-10i).
Overall, these results show that the stochastic a-b ensembles produce a relatively small degree of variability in accumulated precipitation and rain rates. This variability is typically only slightly larger than or comparable to that produced by the θ-pert ensemble but much smaller than from the ICBC or fixed-parameter ensembles. Additionally, there is no clear pattern regarding the chosen autocorrelation scale and its effect on precipitation structure. While STOCH-AB-MID and STOCH-AB-LONG generally produce more spread in domain-accumulated precipitation compared to STOCH-AB-SHORT, this is not the case for rain rate distributions, in which STOCH-AB-SHORT sometimes produces larger variability than the longer autocorrelation scale ensembles.

Cirrus Anvil Properties 5.2.1. Fixed-Parameter Ensembles Compared to ICBC and θ-pert Ensembles
Recall that the stochastic a-b scheme is applied only to unrimed and partially rimed ice, both of which are commonly found in the anvil region of MCSs. Given that the m-D relationship influences size distributions and particle densities, and hence radiative transfer, we explore anvil cloud properties next. Representative plan views of cloud top temperature (CTT) are shown in Figure 11 for SatCORPS (Figures 11a-11c) and WRF CNTL simulations (Figures 11d-11f). The anvil shield is shown at 1300 UTC on 20 May (Figures 11a  and 11d) for the squall line, at 0000 UTC on 24 May after the initiation of supercells (Figures 11b and  11e), and at 0000 UTC on 25 May (Figures 11c and 11f) after the comma cloud associated with the midlatitude cyclone moved into the analysis domain and convection associated with the squall line initiated. Despite differences in location and expanse, all simulations produce an anvil cirrus shield associated with deep convection.
For quantifiable comparison, anvil cloud from observations is identified by locating grid points with CTT < −30°C and a NEXRAD radar reflectivity at 2.5 km above ground level <5 dBZ. Simulated anvil cloud uses the same CTT criteria but requires that the instantaneous rain rate at the lowest model level be <0.1 mm/hr. Simulated rain rates are used instead of reflectivity due to well-established biases in simulated reflectivity (e.g., Stanford et al., 2017;Varble et al., 2011Varble et al., , 2014. SatCORPS optical depth (τ) uses a visible irradiance channel (0.65 μm) and simulated τ is from the Rapid Radiative Transfer Model for General circulation models radiation scheme using SW irradiance between 0.65 and 0.78 μm. Observed ice water path (IWP) is derived from τ and effective particle size, while simulated IWP is computed by vertically integrating ice mass condensate. Since τ uses a visible irradiance channel, simulation-observation comparisons are only performed during time periods in which the domain-mean clear-sky downwelling SW irradiance is larger than 200 W/m 2 (i.e., between 1300 and 0000 UTC). While this misses some of the peak precipitating periods of the 23-25 May cases (see Figure 5), there is still anvil cloud detected during the initiation of supercells on 23-24 May (Figures 11b and 11e) and during the initiation of the squall line and propagation of the comma cloud into the domain during 24-25 May (Figures 11c and 11f).

Journal of Advances in Modeling Earth Systems
All simulations produce greater anvil τ compared to observations for all time periods. This is partly due to the simulations producing larger IWPs compared to observations (not shown), which is likely related to lesser precipitation production (Figure 4). With these biases in mind, some physical inferences can be made by comparing CDFs of anvil τ for the ICBC, θ-pert, FIXED-AB, and FIXED-AB-SAMPLED ensembles. This is shown for three time periods in Figure 12: (1) 1300-2100 UTC on 20 May (Figures 12a-12d), (2) 1300-0000 UTC on 23-24 May (Figures 12e-12h), and (3) 1300-0000 UTC on 24-25 May (Figures 12i-12l).
There is considerable τ distribution variability for the FIXED-AB and FIXED-AB-SAMPLED ensembles. For the 20 May case, FIXED-AB produces median τ spread of up to 55 units (Figure 11c). For the 24-25 May (midlatitude cyclone/squall line) case, median τ spread for FIXED-AB is 70 units (Figure 12k). For all time periods, FIXED-AB generally produces more τ distribution variability than ICBC, and θ-pert produces the smallest variability (Figures 12b, 12f, and 12j). The smaller τ variability for the supercell/bow echo time period on 23 May (Figures 12e-12h) compared to the other two time periods likely results from this case producing less vapor deposition growth, which would be more prevalent in the large anvil regions of the 20 and 24-25 May events that have stronger synoptic forcing associated with upper-level troughs.
We focus on the 20 May time period (specifically, the FIXED-AB ensemble in Figure 12c) to physically interpret the impact of the m-D relationship on τ. The two FIXED-AB members with larger τ CDF values compared to CNTL (yellow and green lines) both have b = 2.7, which acts to shift the mass distribution to

10.1029/2019MS001730
Journal of Advances in Modeling Earth Systems larger particle sizes for a given ice water content (Figure 5a). This mass shift leads to greater size sorting and greater precipitation efficiency, such that after sedimentation of precipitation-sized particles, the ice particles left in the anvil have smaller ice effective radii (r e ) compared to the b = 1.5 simulations. The b = 2.7 simulations also have slower mass-weighted ice fallspeeds (V i ) and larger total ice number concentrations (N i ). Vertical profiles of the FIXED-AB median and mean anvil r e , N i , V i , and total ice mass mixing ratio (q i ) are shown in Figure S1 in the supporting information for the interested reader.
In contrast, the b = 1.5 simulations tend to have geometrically thicker anvils ( Figure S1e) and larger IWP (not shown) that consist of larger particles and smaller N i due to less size sorting and weaker precipitation efficiency. Thus, although the b = 1.5 simulations have thicker anvils, they consist of a relatively small number of larger particles that are not as efficient at extinguishing solar radiation compared to the smaller particles in the b = 2.7 anvils. This reasoning is further supported by the b = 2.7 simulations producing larger τ for the same IWP compared to the b = 1.5 simulations ( Figure S2). Similar trends between anvil optical depth and precipitation following the b-ρ 500μm relationships can be seen in the FIXED-AB-SAMPLED ensemble but are not discussed here.

Stochastic a-b Ensembles and Dependence on Spatiotemporal Autocorrelation Scale
There is clearly a considerable degree of variability in τ distributions produced by the fixed-parameter ensembles shown in Figure 12. While the stochastic ensembles are not expected to produce the same degree of variability and the same physical inference cannot be drawn due simply to the stochastic nature, Figure 13 shows that the stochastic scheme does alter the τ distribution. For example, STOCH-AB-MID produces median τ variability of 15 units for the 23-24 May case ( Figure 13e) and of 35 units for the 24-25 May case (Figure 13h). STOCH-AB-LONG for both cases has a spread in median τ of 10-15 units.
Variable τ distributions among stochastic ensemble members can occur through variability in IWP or r e passed to the radiative transfer scheme. Average τ values as a function of IWP for the stochastic a-b ensembles and the ICBC ensembles are shown in Figure 14. The relationship between average τ and IWP is nearly constant among all ICBC members, regardless of the simulated case, and does not vary from case to case (Figures 14d, 14h, and 14l), consistent with almost no spread in profiles of r e (not shown). In contrast, there is considerable variability in the τ-IWP relationship among the stochastic a-b members, because the m-D

Journal of Advances in Modeling Earth Systems
relationship is perturbed which directly influences r e . For example, two members of STOCH-AB-LONG for the 20 May case (blue and orange lines) produce higher τ than observed for a given IWP, while three members produce less, enveloping the observed relationship ( Figure 14c). These two members that produce larger τ also produce smaller r e profiles, slower V i , and larger averaged b parameters compared to the other three simulations, which is in agreement with the FIXED-AB ensemble showing larger τ for a larger b parameter. In the STOCH-AB-MID ensemble for the 24-25 May time period, one member produces a relationship very close to what was observed (Figure 14j). Note also that the observed τ-IWP relationship is different for each case and that the stochastic scheme is better able to capture this variability. Overall, these results show that the stochastic a-b scheme produces τ distribution variability for a given IWP, while this variability is virtually nonexistent for the ICBC ensembles.
The τ distribution variability shown in Figure 13 conceivably impacts simulated SW cloud radiative forcing (CRF). We show only the 20 May case for brevity, but results for the other case are qualitatively similar. The time series of mean anvil SW CRF (i.e., mean among all columns identified as anvil) for the 20 May case is shown in Figure 15. Here, CRF is defined as the difference between the simulated SW downwelling irradance and the clear-sky SW downwelling irradiance at the surface. Observations are shown in black, in which the observed mean SW CRF is the average of the radiometer stations shown in Figure 2 at times when the stations were under anvil cloud coverage. Note that several thousand anvil columns are present in the simulation means for a given time, while the observations include at most 13 values. Observed standard errors around the mean are included in Figure 15. While the small number of radiometers limit quantitative comparison, some inferences can be made.
All simulations produce greater anvil SW CRF compared to observations (Figure 15), which is in agreement with simulations producing higher τ compared to observations ( Figure 13). Nonetheless, simulations show SW CRF spread of up to 50 W/m 2 at 1800 UTC for the STOCH-AB-LONG ensemble (Figure 15c). In

Journal of Advances in Modeling Earth Systems
addition, the two simulations that produce the highest τ in the STOCH-AB-LONG ensemble (blue and orange lines in Figures 13c and 14c) also produce the largest SW CRF in Figure 13c and have the smallest ice effective radii (not shown). As expected, there is a link between CRF, τ, and ice effective radiusensemble members with larger CRF are linked to larger τ and smaller ice effective radius profiles. Similar consistency between τ distributions, anvil SW CRF, and ice effective radius are seen in the STOCH-AB-SHORT and STOCH-AB-MID ensembles.
Anvil SW CRF normalized by the clear-sky downwelling SW irradiance is shown in Figures 15d-15f, which removes the diurnal cycle. These results indicate that anvil SW CRF varies for the stochastic simulations by a few percent across the entire time period, reaching up to 5% at some times for the STOCH-AB-LONG ensemble (Figure 15f).

Stochastic Ice Fallspeed Ensembles 5.3.1. Precipitation Structure
Additional simulations are run in which the ice fallspeed (V i ) is stochastically perturbed, as described in section 2.3. Motivation for this sensitivity test results from limited variability in precipitation properties using the stochastic m-D scheme. We note that altering the m-D relationship of unrimed and partially rimed ice alters V i since the fallspeed is proportional to the ratio of mass and cross-sectional area. However, the stochastic V i simulations are applied to all ice (including rimed ice) since this is expected to more significantly impact precipitation structure in convective system simulations. Because rimed ice as well as unrimed and partially rimed ice are perturbed in the stochastic V i scheme, sensitivity to this scheme should not be compared directly to that from the stochastic m-Dscheme.
Ensemble mean domain-accumulated volumetric precipitation for the fixed V i and stochastic V i ensembles show that for both cases there is no significant shift in mean behavior of precipitation evolution (Figures 16a and 16d). However, there is variability produced by the FIXED-VI and FIXED-VI-SAMPLED ensembles that at times exceeds the variability produced by ICBC (Figures 16b, 16c, 16e, and 16f).
Although not shown, patterns of accumulated precipitation produced by the FIXED-VI and FIXED-VI-SAMPLED ensembles behave as expected as a function of the multiplicative factor, in which larger fallspeeds lead to more accumulated precipitation and smaller fallspeeds lead to less (which is the reason for larger σ shown in Figures 16c and 16f). Interestingly, the FIXED-VI member with the multiplicative factor of 0.5 produces the least accumulated precipitation and smallest areas of moderate-to-high hourly rain rates but the largest total precipitating area. On the other hand, the FIXED-VI member with the largest multiplicative factor of 2 produces the most accumulated precipitation, larger areas of moderate-to-high hourly rain rates, but the overall smallest total precipitating area. Thus, increasing the fallspeed increases precipitation maxima but decreases the total area of precipitation.
The stochastic V i scheme is able to produce larger σ compared to the θ-pert ensemble for all time periods (Figures 16b,16c,16e,and 16f). This is not surprising given studies that indicated the large impact of ice fallspeed on deep convection (e.g., Bryan & Morrison, 2012;McCumber et al., 1991;Morrison & Milbrandt, 2011). While the stochastic V i scheme is not expected to produce the same variability as the fixed-parameter

Journal of Advances in Modeling Earth Systems
V i ensembles, these results suggest that allowing for variable V i can alter condensate fluxes that impact surface precipitation. Additionally, the 20 May case exhibits some dependence on spatiotemporal autocorrelation scale, in which STOCH-VI-SHORT produces the smallest variability of the stochastic V i ensembles while STOCH-VI-LONG and STOCH-VI-MID produce larger σ. This is true for the 23-25 May case as well. Also note that normalized σ in the 23-25 May case is much larger compared to the 20 May case, with normalized σ in the stochastic ensembles frequently exceeding 10%.
Rain rate contributions to total rainfall CDFs for the stochastic V i ensembles and the FIXED-VI ensembles are shown in Figure 17. FIXED-VI produces considerable variability (Figures 17d, 17h, and 17l), consistent with the expected relationship of increasing V i leading to higher rain rates. There are several instances where the stochastic V i ensembles contain members in which the rain rate distribution is closer to observed. For example, STOCH-VI-MID during the 23-24 May supercell/bow echo event produces members that generally envelop the ABRFC distribution ( Figure 17f). All stochastic V i ensembles during the 24-25 May midlatitude cyclone/squall line case produce median rain rate variability of 1-2 mm/hr, with the stochastic fallspeed ensembles generally producing rain rates closer to observed compared to CNTL (Figures 17i-17k). The FIXED-VI ensemble also produces variability comparable to or greater than the ICBC ensembles (cf. Figure 17 and Figure 9).

Cirrus Anvil Properties
The stochastic V i scheme also produces variability in τ distributions but for different reasons compared to the stochastic m-D scheme. CDFs of IWP and τ are shown in Figures 18a and 18b, respectively, for the 20 May STOCH-VI-LONG ensemble. Clearly, the same members producing higher τ also produce higher IWP, and this tight τ-IWP relationship is seen in Figure 18c. Members with higher τ and IWP have higher Thus, CRF variability for the stochastic fallspeed ensembles largely results from variability in IWP. This is also implied by a lack of variability in ice effective radius profiles (not shown). We note that τ variability for the stochastic m-D ensembles ( Figure 13) also partially results from variable IWP distributions, but IWP variability for the stochastic a-b ensembles is less than in the stochastic V i ensembles. However, as discussed earlier, τ variability for the stochastic m-D simulations exists even for the same IWP ( Figure 14) because of variability in ice effective radius profiles (the reasoning for which is discussed above). Overall, these results suggest that CRF variability for the stochastic m-D ensembles results from spread in both effective radius and IWP, while CRF variability for the stochastic fallspeed ensembles is mostly from IWP spread alone.

Discussion and Conclusions
A stochastic framework was implemented into the P3 microphysics scheme (Milbrandt & Morrison, 2016;Morrison & Milbrandt, 2015), in which parameters of the unrimed and partially rimed ice mass-size (m-D) relation vary spatially and temporally during simulations. This addresses a current lack of knowledge regarding the impacts of natural variability of ice particle properties on cloud and precipitation properties. Here, we varied the m-D relationship, which has an established observational variability (e.g., Finlon et al., 2018). In addition, we explored the sensitivity to parameter uncertainty via ensembles with parameters that are fixed in time and space but varied among simulations (so-called "fixed-parameter" ensembles). We note that the stochastic m-D framework was applied only to unrimed and partially rimed ice and not to graupel or small, spherical ice.
Although a few studies have investigated the impact of including stochastic microphysics in convection simulations (e.g., Kuell & Bott, 2014;Qiao et al., 2018), our approach is unique in several ways. First, by guiding the sampling distribution of stochastic parameters from aircraft observations, parameter values remained within physical limits of what is observed in nature. Second, we compared the stochastic scheme against a suite of other ensemble configurations, which was not previously investigated in the context of stochastic microphysics parameterization.
Observational guidance for the scheme was provided by aircraft measurements obtained during the MC3E field campaign. The stochastic scheme was supplied with m-D parameter distributions retrieved from these observations and sampled according to a prescribed spatiotemporal autocorrelation scale. Simulations were performed for two deep convection cases. In addition to a control run employing the nonstochastic P3 scheme, the stochastic simulations used various spatiotemporal autocorrelation scales and were compared against fixed-parameter ensembles, ensembles employing large-scale initial and lateral boundary condition perturbations (ICBC), and an ensemble with only small-amplitude "noise" added to the lower-level potential temperature field.
Domain-accumulated volumetric precipitation was only slightly affected by the stochastic m-D scheme, with slightly more spread compared to the grid-scale noise ensemble but with no shift in ensemble mean behavior

Journal of Advances in Modeling Earth Systems
closer to the observed amount. More spread relative to ensemble mean was produced during the bow echo event on 24 May compared to the 20 May squall line case, which is expected given the 20 May case was strongly forced at the synoptic scale, whereas the bow echo was not. Compared to the stochastic ensembles, larger variability was produced by the fixed-parameter ensembles, indicating a sensitivity to fixed-parameter values that is much larger than the sensitivity to parameter variability. The stochastic scheme also produced some rain rate distribution variability. Median rain rates in the stochastic ensembles varied by up to 1.5-2 mm/hr, though at other times the variability was indistinguishable from the grid-scale noise ensemble. Modest precipitation structure variability produced by the stochastic m-D scheme and greater variability produced by the fixed-parameter ensembles were both much smaller than variability produced by the ICBC ensemble.
Analysis of cirrus anvil properties showed that the stochastic m-D scheme had important impacts on radiative transfer. The fixed-parameter ensemble in which the m-D parameters were chosen at the extremes of the observationally guided range had variability in median optical depth (τ) up to 70 units. The stochastic m-D ensembles produced less spread in the τ distributions than the fixed-parameter ensembles, which is not surprising but still had a spread in median τ of up to 15 units. Stochastic ensemble members that produced larger (smaller) τ also produced profiles of smaller (larger) mean effective radius and greater (lesser) cloud radiative forcing (CRF). There was variability in τ for a given IWP, so that the stochastic m-D ensembles had variability in radiative transfer (thus altering CRF) even for similar vertically integrated anvil ice mass profiles. In contrast, the relationship between τ and IWP did not vary for the ICBC ensemble, even though it produced the largest spread in τ distributions. Furthermore, the observed relationship between τ and IWP varied from case to case. The stochastic m-D scheme was able to capture this variability, while the nonstochastic simulations could not.
Comparison of the fixed m-D and V i ensembles with the stochastic ensembles can be interpreted from the standpoint of systematic versus random model error. Both the fixed-parameter and stochastic ensembles were sampled from the same parameter distributions, with the implicit assumption that the amplitude of the systematic and random errors are equal (note that the fixed-parameter ensembles are analogous to the stochastic ensembles with a parameter perturbation spatiotemporal autocorrelation scale much larger than the model domain size and integration time). Results from these ensembles suggest that cloud and precipitation structure is more sensitive to systematic model error (via the fixed-parameter ensembles) than random error (via the stochastic ensembles). This is not surprising as systematically changing the m-D and parameters in the fixed-parameter ensemble forces the unrimed and partially rimed ice particles to take a specific value of bulk particle density, resulting in systematic differences in model output, rather than allowing the density to vary in time and space as in the stochastic ensembles. A similar situation holds for perturbing V i .
Additional ensembles were analyzed in which the prefactor of the ice fallspeed-size relation (V=cD d ) was varied by a multiplicative factor with ±2 standard deviations from the mean, lying between 0.5 and 2. Results showed that the stochastic V-Dscheme produced larger variability in rain rate distributions compared to the stochastic m-D scheme, often producing higher rain rates closer to what was observed. The stochastic V-D scheme also produced spread in τ and CRF distributions, but the τ spread was highly correlated with the IWP spread, and virtually no variability was produced in the τ-IWP relationship. Thus, τ and CRF variability produced by the stochastic V-D scheme resulted almost entirely from differences in IWP owing to impacts on sedimentation. In contrast, τ and CRF variability produced by the stochastic m-D scheme was caused by spread in both IWP and ice effective radius.
Currently, there is little observational guidance of appropriate spatiotemporal autocorrelation scales for the parameters that were tested. Next steps should focus on deriving these scales in a physically based manner from observations. Nonetheless, we examined scales ranging from individual clouds O(10 min, 10 km) to effectively an infinite spatiotemporal scale (i.e., the m-D parameters were randomly sampled from the observationally guided distribution but held fixed in time and space). The shortest spatiotemporal scale generally produced the smallest anvil cirrus property variability (similar to that of the grid-scale noise ensemble), whereas the sampled, fixed-parameter ensemble produced the most variability. This suggests that increasing the spatiotemporal autocorrelation scale increases precipitation structure and anvil cloud property variability. However, more research is needed to determine appropriate autocorrelation scales from observations and how they relate to the results presented here.

Journal of Advances in Modeling Earth Systems
We point out several limitations of this study. In particular, systematic errors leading to underpredicted precipitation rate are likely associated with initial and lateral boundary condition uncertainty not captured within the ICBC ensemble. However, these systematic errors are unlikely to affect the conclusions drawn here, as results were broadly consistent among the multiple events simulated. The scheme also did not include explicit vertical variability or the temperature dependence established in Finlon et al. (2018), which may be potentially important to include in future work. In addition, the stochastic framework has only been implemented into a single microphysics scheme and for only a few parameters. Future work with stochastic microphysics should explore additional parameters, be extended to other microphysics schemes, and tested for additional cases of varying convective mode and thermodynamic environments. For instance, the stochastic m-D scheme may produce more variability in simulated weather systems in which the precipitation is dependent more on diffusional growth rather than riming. A stochastic framework that fully accounts for variable particle shape should also include variability of the area-dimension relationship that retains consistency with the the m-D parameters for the fallspeed calculations. Finally, it remains unclear whether the stochastic microphysics framework will increase skill in numerical weather prediction models. In order to test this, simulations should be run on several forecast cycles (i.e., over a season) and compared with other ensemble configurations (e.g., multiphysics, SPPT, stochastic kinetic energy backscatter, and ICBC perturbations).