Sensitivity of Modeled Microphysics to Stochastically Perturbed Parameters

This study examines the characteristics of several model parameter perturbation methodologies for ensemble simulations of cloud microphysical processes in convection. A simplified 1D model is used to focus the results on cloud microphysics without the complication of feedbacks to the dynamics and environment. Several parameter perturbation methods are tested, including non‐stochastic and stochastic with various distributions and parameter covariance. We find that an ensemble comprised of different time‐invariant parameters (non‐stochastic) exhibits little bias, but small spread. In addition, its behavior does not respect the time evolution of convection through its various phases. Stochastic parameter (SP) methods in which no inter‐parameter covariance is applied produce greater spread, but significant bias. The bias is particularly large for lognormal parameter perturbation distributions. The ensemble spread is retained and the bias reduced when time‐varying parameter covariance is applied. In this case, the SP scheme is able to adapt to the time and state‐dependent covariance structures and produce ensemble characteristics that are consistent with the specific microphysical processes operating at any given time. The results suggest that SP schemes would benefit from inclusion of parameter covariances, and specifically those that vary with the state of the system. It also suggests that a Normal or LogNormal SP scheme with no covariance may significantly impact the ensemble bias. Finally, the results indicate that high temporal and spatial resolution observations may be needed to characterize the variability in parameter values and covariance.

Over the last two decades, the inclusion of stochastic representation of unresolved variability of sub-grid scale processes in weather and ES models has emerged as a viable approach to representing the intrinsic uncertainties of physical parameterizations in the context of improving probabilistic prediction skill using ensembles (Berner et al., 2011(Berner et al., , 2015Bouttier et al., 2012;Buizza et al., 1999;Charron et al., 2010;Christensen et al., 2015;Davini et al., 2017;Hermoso et al., 2021;Jankov et al., 2017;Leutbecher et al., 2017;Palmer et al., 2005Palmer et al., , 2009Ollinaho et al., 2017;Romine et al., 2014;Sanchez et al., 2015;Stanford et al., 2019;Thompson et al., 2021). The representations of uncertainty using this approach are commonly referred to as stochastic parameterizations . Stochastic parameterizations may be formulated either by including an empirical stochastic representation of the parameterization uncertainty into a numerical model a posteriori Palmer, 2001), or by derivation of stochastic governing equations (Debussche et al., 2012). The latter formulation Abstract This study examines the characteristics of several model parameter perturbation methodologies for ensemble simulations of cloud microphysical processes in convection. A simplified 1D model is used to focus the results on cloud microphysics without the complication of feedbacks to the dynamics and environment. Several parameter perturbation methods are tested, including non-stochastic and stochastic with various distributions and parameter covariance. We find that an ensemble comprised of different time-invariant parameters (non-stochastic) exhibits little bias, but small spread. In addition, its behavior does not respect the time evolution of convection through its various phases. Stochastic parameter (SP) methods in which no inter-parameter covariance is applied produce greater spread, but significant bias. The bias is particularly large for lognormal parameter perturbation distributions. The ensemble spread is retained and the bias reduced when time-varying parameter covariance is applied. In this case, the SP scheme is able to adapt to the time and state-dependent covariance structures and produce ensemble characteristics that are consistent with the specific microphysical processes operating at any given time. The results suggest that SP schemes would benefit from inclusion of parameter covariances, and specifically those that vary with the state of the system. It also suggests that a Normal or LogNormal SP scheme with no covariance may significantly impact the ensemble bias. Finally, the results indicate that high temporal and spatial resolution observations may be needed to characterize the variability in parameter values and covariance.
Plain Language Summary All numerical weather and climate prediction models contain approximate representations (parameterizations) of high resolution cloud processes. Model error results from either inaccurate parameter settings, or from the fact that parameters are static while they should vary in time and space. Previous studies have shown that random variation of model parameters can help to improve predictability of weather and climate. This study examines several methods by which model parameters can be varied in time, and finds that the best results are obtained when parameter variations are consistent with the time evolution of the weather system that is being simulated.
The most commonly used empirical stochastic parameterization, that is directly related to the uncertainty of the existing physical parameterizations in numerical models, is the Stochastically Perturbed Parameterization Tendency (SPPT; Palmer et al., 2009) scheme. The scheme consists of adding random perturbations to net thermodynamic and dynamic tendencies produced by physical parameterizations at resolvable scales. The perturbations are sampled at every time step by multiplying the unperturbed net tendency by a stochastic temporally and spatially correlated noise. Spatio-temporal variability of the noise is governed by a linear Auto Regressive process model (AR1) applied to a Gaussian correlated spatial noise pattern in the horizontal domain. Applications of the SPPT scheme in global models have produced significant positive impacts on ensemble prediction skill by increasing the ensemble spread and reliability and reducing biases in some variables for prediction ranges from days to seasonal climate (Bouttier et al., 2012;Charron et al., 2010;Leutbecher et al., 2017;Ollinaho et al., 2017;Palmer et al., 2009;Sanchez et al., 2015;Weisheimer et al., 2014). The most significant impacts were found in the tropical regions associated with convective and precipitation processes . Additionally, positive impacts have been demonstrated in climate model simulations for processes such as the Madden-Julian Oscillation , El Nino Southern Oscillation , and tropical precipitation (Davini et al., 2017). Studies with regional high-resolution models using the SPPT scheme have shown improvements also to short-range ensemble prediction, most notably to the ensemble spread (Baker et al., 2014;Berner et al., 2011Berner et al., , 2015Jankov et al., 2017Jankov et al., , 2019Romine et al., 2014;Wang et al., 2019).
Although including the SPPT scheme as a stochastic representation of parameterized physics uncertainty in deterministic numerical models has shown clear benefits, it has also been noted that the SPPT could lead to numerical instability, and to physical inconsistencies and imbalances in energy and moisture fluxes locally and globally (Davini et al., 2017;Leutbecher et al., 2017;Ollinaho et al., 2017). These shortcomings have been attributed to the assumption of a uniform uncertainty representation for all parameterizations in the model, implying the same error characteristics for different processes, and to the absence of conservation constraints on the tendency perturbations . To address the need for a more physically-based approach, recent developments have focused on either the implementation of stochastic schemes for specific processes, such as for example, the convective initiation in boundary layer (Kober & Craig, 2016) and deep convection (Bengtsson et al., 2011(Bengtsson et al., , 2019, or on holistic formulations that may apply to different parameterizations by using a single formulation for stochastic perturbation of selected parameters within the deterministic parameterizations (Baker et al., 2014;Christensen et al., 2015;Jankov et al., 2017;McCabe et al., 2016;Ollinaho et al., 2017;Thompson et al., 2021). Within the latter framework, hereafter referred to as the Stochastic Parameters (SP) approach, a parameter is considered to be any tunable scalar physical quantity within an existing deterministic parameterization. In the current study we consider the SP approach, focusing on cloud and precipitation microphysics parameterizations.
The mathematical formulation of an SP scheme is similar in simplicity to an SPPT scheme, in that they both make use of a single stochastic process model for any parameterization, but SP applies it to the parameters within each selected parameterization instead of to the tendencies. In this way, the sources of uncertainty unique to each parameterization are represented while preserving the parameterization's built-in physical consistency (Baker et al., 2014;Ollinaho et al., 2017). As a result, the SP schemes produce physically-based and physically-consistent stochastic variability to the tendencies of the resolved states. Evaluations of the SP schemes in global and regional models have shown positive impacts on the ensemble prediction skill of various degrees depending on the prediction time range, the modeling system, and verification variables (Christensen et al., 2015;Hermoso et al., 2021;Jankov et al., 2017;Leutbecher et al., 2017;McCabe et al., 2016;Ollinaho et al., 2017;Stanford et al., 2019;Thompson et al., 2021;Wang et al., 2019). In comparison to the SPPT, the impacts of the SP schemes on ensemble spread were found to be smaller, but adding value when the two approaches were combined (Christensen et al., 2015;Hermoso et al., 2021;Jankov et al., 2019;Ollinaho et al., 2017). These findings point to two non-exclusive conclusions: the unresolved variability has a broader scope than is contained in the parameter uncertainty, and applications of the SP approach may need to capture a larger diversity of processes (i.e., more parameters and parameterizations within a model) with improved process-specific configurations and better understanding of impact mechanisms to achieve full benefits (Jankov et al., 2019;Leutbecher et al., 2017). The Recent studies on the application of the SP schemes to microphysics parameterizations using cloud-permitting resolution regional models provided examples of benefits from including the uncertainty of these processes. Stanford et al. (2019) showed that perturbing the ice microphysics parameters involved in mass-size and fallspeed-size relationships using the SPP (Stochastically Perturbed Parameterization; Ollinaho et al., 2017) scheme resulted in significant impacts on the variability of ice cloud optical depth, rain-rate and precipitation in ensemble simulations of deep convection events. Similarly, Hermoso et al. (2021) demonstrated that the Random Parameters (RP) scheme (McCabe et al., 2016) applied to several microphysics parameters in a two-moment bulk microphysics scheme significantly increased ensemble variability of convective storm evolution for the studied case. Thompson et al. (2021) reported important localized impacts on the precipitation properties including changes to maximum hail-size from perturbing only two microphysics parameters using the SPP scheme. These results clearly point to potential for expanding the scope of relevant impacts by including more processes within the SP framework. Another direction for further development that so far has not received much attention concerns configurations of the SP schemes that involve selection of controlling factors, such as temporal and spatial decorrelation scales and statistical distributions for the parameters. The evaluation and calibration of the SP configurations in the prior studies have focused mostly on enhancing the ensemble prediction skill, in particular the ensemble spread (Christensen et al., 2015;Jankov et al., 2019;Ollinaho et al., 2017;Thompson et al., 2021). Although justified for the purpose, this approach does not consider representativeness of the configurations with respect to the process-level uncertainty that the schemes are meant to capture. This could lead to selections of the controlling factors that are inconsistent with the properties of the physical processes represented by the parameterizations .
This study addresses evaluation of the SP configurations at the process level, focusing on cloud microphysics parameterizations. Our approach is based on the notion that the parameterization uncertainty is a statistical generalization of the parameterization errors with respect to the true processes and as such could be characterized through evaluation against reference data with a high information content about the processes. Candidate approaches include: using Cloud Resolving Models (CRMs) including Large Eddy Simulation configurations, with more detailed representation of cloud microphysics, and comparison versus observations with high information content about process level outputs from the microphysics parameterizations. Both approaches have strengths and weaknesses. The CRMs produce dynamically consistent high information content data but carry their inherent numerical model uncertainties and due to computational burden are also challenging to apply to a wide range of natural scenarios. They have been used for evaluation and calibration of microphysics parameterizations (Khain et al., 2015) and may be in principle adopted to characterize properties of SP-versions of the parameterizations. That methodology would be similar to the coarse-graining method that was applied to evaluate and propose modifications to the SPPT scheme using CRM simulations (Shutts & Pallares, 2014). We are not aware, however, of any published studies that have used CRM data to calibrate SP-microphysics schemes. The comparison versus observations approach has the strength of using data that represent the actual true state with high accuracy (to the limits of observation errors) and for many scenarios. The primary weakness of this approach is that the data are by design limited to the observable quantities, requiring careful selection of types of observations that could provide sufficient information about the process-level state of the microphysics. In this study the evaluation of the SP-microphysics configurations was considered from the prospective of observability by satellite remote sensing.
As a first step we investigated sensitivity of a microphysics parameterization to different configurations of the SP scheme applied to multiple microphysical parameters within a Lagrangian column model that was used in prior studies to characterize parameter uncertainty by means of multivariate nonlinear inversions using simulated remote sensing observations (Posselt & Vukicevic, 2010;van Lier-Walqui et al., 2012. The column microphysics model is forced with prescribed time-varying profiles of temperature, humidity and vertical velocity. This modeling framework allowed for investigation of the effect of changes in microphysics parameters on the model output in isolation from any feedback to the cloud-scale dynamics. The test case selected in this study of an idealized representation of mid-latitude squall-line convection is the same as was used in prior studies (Posselt & Bishop, 2012;Posselt et al., 2014;Posselt & Vukicevic, 2010;van Lier-Walqui et al., 2012. Use of this model in prior studies enabled use of multi-parameter distributions from previous Bayesian inversion as the basis for setting the statistical properties of parameter distributions used in the SP scheme. Additionally, impacts of the non-stochastic and stochastic multi-parameter representation of parameterization uncertainty on the microphysics model solution and the related satellite-based observable quantities could be directly compared.
The current study addresses the following questions.
What is the microphysics state response to stochastically perturbed parameters ?
What is the sensitivity of observable quantities to the SP configurations ?
What is the impact of using correlated parameter perturbations ?
What is the potential of using microphysics sensitive observations to estimate the SP configuration controls ?
The paper is organized as follows. The 1D column model, the formulation of the SP scheme for the study of configurations, and the design of the sensitivity experiments is presented in Section 2. The experiment results are discussed in Section 3. The summary and conclusions are outlined in Section 4.

1D Column Microphysics Model and Case
Our goal is to explore stochastic parameter perturbations in a system that is low-dimensional and simple, enabling straightforward interpretation of results, and also sophisticated enough to represent the complexity of a realistic three dimensional cloud system. Posselt and Vukicevic (2010) described a single column version of the NASA Goddard Cumulus Ensemble (GCE) model (Tao et al., 2014) that simulates the changes a vertical column of the atmosphere would encounter if it were to be embedded in a three dimensional deep convective system. By applying a base state temperature and water vapor profile representative of a deep convective environment, then imposing perturbations to temperature, water vapor, and vertical motion consistent with convective system evolution from initiation through maturity to dissipation, the model can be made to approximate the behavior of a single column within a fully three dimensional environment. The model microphysical parameterization is consistent with those used in modern numerical weather prediction models, and is a single moment bulk scheme (Lang et al., 2007;Lin et al., 1983;Rutledge & Hobbs, 1983, 1984Tao et al., 2019) with two liquid (cloud droplets and rain) and three ice (suspended ice particles, falling snow/aggregates, and graupel/hail) hydrometeor species. Clouds interact with longwave and shortwave radiation, and are allowed to be advected in the vertical direction only. No sources or sinks of cloud mass due to horizontal advection are considered.
In the setup described in Posselt and Vukicevic (2010), the model emulates the development of a linear squall line type of convection over a three-hour time period, represented by time evolution of the vertical profiles of cloud droplets, rain, pristine ice, graupel and snow as shown in Figure 1 (respectively panels g-k). This is desirable for our experiments, as the evolution proceeds through two distinct phases: convective (0-90 min), during which there is more intense vertical motion and also more intense rainfall, and stratiform (120-180 min), during which there is more moderate vertical motion and lighter rainfall. The 90-120 min time period represents a transition period in between the two phases. In addition, different cloud microphysical processes are active in convective versus stratiform time periods, with convective times controlled more strongly by liquid microphysics and stratiform times controlled more by ice microphysical parameters. This change in behavior allows us to test how the SP scheme is able to adapt to changes in the degree to which various parameters influence the model output.
The specific parameters we perturb (Table 1) control: the snow and graupel particle fall speed, snow and graupel density, rain particle size distribution, and cloud-to-rain autoconversion threshold. As in Posselt and Vukicevic (2010), we use as output quantities of interest the Precipitation Rate (PR), Liquid and Ice Water Paths (respectively LWP and IWP), and Outgoing Longwave and Shortwave (respectively OLW and OSW) radiation, all output from the model as snapshots at regular times. The model outputs shown in Figure 1 were used as the reference control simulation in this study.

Sampling of Stochastically Perturbed Parameters
The commonly used SP schemes for the global and regional ensemble models are based on the linear AR1 models. They differ however in the selection of the AR1 model class (Grunwald et al., 2000), the quantity to which the model is applied and the underlying statistical parameter distribution used for sampling the innovations at every time step. For example, the SPP scheme first formulated by Ollinaho et al. (2017) and since applied in a number of other studies makes use of the Gaussian AR1 model applied to spectral coefficients of 2D horizontal Gaussian correlated noise with the innovations based on either the Log-Normal or Normal univariate distributions (Jankov et al., 2017(Jankov et al., , 2019Ollinaho et al., 2017;Stanford et al., 2019;Thompson et al., 2021). Another widely used SP scheme referred to as the Random Parameters (RP) uses a non-Gaussian linear AR1 model applied directly to spatially invariant parameters, with the innovation term based on the bounded uniform distribution (Hermoso et al., 2021;McCabe et al., 2016). In this study the generalized formulation of the Gaussian linear AR1 model is used (Grunwald et al., 2000) from which these commonly used schemes could be derived as special cases as shown in the Appendix A.
The Gaussian linear AR1 model is expressed Here Φ is a vector stochastic quantity with mean Φ , is the first order regression coefficient determined by the model time step and decorrelation time scale Λ , 1∕2 Φ ( ) is the innovation term where 1∕2 Φ is square root covariance matrix, and ( ) is a (0, ) random vector, with being the identity matrix. This model constitutes the SP scheme whose configurations would depend on the underlying sampling distribution for the innovations, the associated first and second moment statistics, and the decorrelation time scale Λ .
The Log Normal (LN) and Normal (N) distributions are considered for representing the variability of microphysics parameters within a range of permissible values (Table 1) , where ′ and are respectively the perturbed and unperturbed parameter values, the latter being the constant parameter value used in the deterministic microphysics model. Assuming for simplicity that the parameters are uncorrelated would make a random ( 2 ) quantity. Using the condition that the mean of the perturbed parameter is equal to the unperturbed value renders = − 2 2 and 2 = , where 2 is the assumed (i.e., the prescribed) parameter variance. Using the N distribution would give = ′ with the mean = and the variance 2 . When the parameters are assumed correlated a non-diagonal 1∕2 Φ would need to be specified for either choice of the distribution. Concerning the influence of the decorrelation time scales in this study the relationship = − Λ is used.
For the sensitivity experiments the second moment statistics for the LN and N parameter distributions were determined based on the data obtained in the prior study, where the uncertainty of microphysics parameters in the 1D column model was estimated by MCMC (Markov Chain Monte Carlo) inversions using simulated satellite-based observable quantities of PR, IWP, LWP, OSW and OLW radiation (Posselt & Vukicevic, 2010). In that study, the observations 30 min apart were derived from the 1D column model solution equal to the one shown in Figure 1 at 30, 60, 90, 120, and 150 min. The MCMC inversions were performed for 10 microphysical parameters, resulting in a 10-dimensional joint empirical statistical distribution.
The data from the prior study (https://doi.org/10.5281/zenodo.5736940, Posselt et al., 2021) consisting of over 1 million samples of valid parameter vectors (the parameter combinations that correspond to the model solution within the assumed error range for the simulated observations) were used to compute variances and square root covariance matrices via the Cholesky decomposition for subsets of the 7 and 8 parameters considered in this study for the sensitivity experiments. Initial sensitivity experiments using either the LN or N distribution showed Note. The standard deviations in the last column correspond to 20% of the values obtained using the inversion data from Posselt and Vukicevic (2010) study.

Table 1 Microphysics Parameters Selected for Perturbations in the Sensitivity Experiments With the Corresponding Prescribed Mean Values, Range Boundaries and Standard Deviations Used for Sampling the Parameter Variations by Different Configurations of the SP Scheme
that the parameter variances based on the inversion data were too high for the application with these distributions in the SP scheme; they resulted in a high probability of perturbed parameter values that were close to the theoretical boundaries of the parameter ranges. For this reason, the variances were reduced to 20% of the original values. The parameters selected for perturbation in this study together with the corresponding units, value ranges, mean and standard derivations are presented in Table 1. The mean values were set to the "true" model solution that was used to simulate the observations in the prior study. The inversion data were used also to randomly select vector parameter samples for ensemble simulations assuming time-invariant parameters. This configuration was used to compare the stochastic versus non-stochastic parameter perturbations. The non-stochastic configuration is equivalent to using the multi-parameter ensemble strategy (Yussouf & Stensrud, 2012).

Sensitivity Experiments
The sensitivity experiments with the SP scheme involved 100-member ensemble simulations, in which each member was evolved with a different stochastic sequence of parameter perturbations. Additionally, the same size ensemble simulations with randomly sampled constant-in-time parameters described in the prior section were performed, referred to as non-SP configuration. The SP experiments explored impacts of using different statistical distributions (LN and N) for sampling the innovations, uncorrelated and correlated parameter perturbations, and different decorrelation time scales. The comparison of correlated versus uncorrelated perturbations was performed using only the N distribution. All configurations were applied to 7 and 8 parameter sets and decorrelation times of 1, 3, 6 and 12 hr (Table 2).
Two sets of sensitivity experiments were performed using the correlated parameter sampling: one with time-invariant and other with time-varying covariance. The time-invariant covariance was computed using the prior study data as discussed in the previous section. The data for the time-varying covariance were generated in this study motivated by the initial sensitivity experiments, which showed significant variability of the degree of model sensitivity to the SP scheme depending on the storm phase. The results from the sensitivity experiments are discussed in the following section. Here, we present the approach taken to estimate the storm-phase-dependent covariances. New MCMC inversions were performed using the same algorithm as in Posselt and Vukicevic (2010) for three different observing periods instead of the five observation instances (30, 60, 90, 120, and 150 min) within one three-hour model integration as was done in the prior study. The observing periods were set based on the phases of the simulated storm evolution identified as: convective (30, 60, 90 min), transition to stratiform (90, 120 min) and stratiform (120, 150, 180 min). The new inversions used the same observable quantities as in the prior study (PR, LWP, IWP, OSW and OLW), but three separate observing time periods. This resulted in three different empirical joint parameter distributions. Each corresponding covariance matrix was then applied sequentially in the SP scheme configuration with the time-varying covariance (labeled SP-NCovV in Table 2).
Comparison of correlation matrices corresponding to the prior and new inversion data ( Figure 2) indicated significant differences in the relationships between the two representations of the covariance parameter variations for certain parameters. For example, the correlation based on the entire storm evolution (the prior study) between Note. PV10 refers to the origin of the inversion data obtained in the prior study (Posselt & Vukicevic, 2010). Experiments with 7 and 8 perturbed parameters included, respectively, the first 7 and all parameters shown in Table 1.

Table 2
Parameter Sampling Configurations Used for the Sensitivity Experiments the coefficients involved in the fall-speed parametrization for the snow and graupel (respectively and ) has significantly higher amplitude (+0.47) and is of the opposite sign than the correlations computed from the new data for either phase of the storm (respectively Figures 2a-2d). Positive correlation between snow and graupel fall speed coefficients results from the constraint that < in the MCMC experiments. As increases, must necessarily increase as well. The correlations corresponding to the different storm phases are all of small amplitude and negative (Figures 2b-2d) with the smallest amplitude during the convective phase (−0.06; Figure 2b) and a gradual increase toward the stratiform phase (Figures 2c and 2d). While the parameter value constraint was also enforced for subset time intervals, the snow and graupel parameters had more limited influence relative to the entire simulation. Also for the parameters involved in the rain processes, the threshold mass mixing ratio for auto-conversion to rain ( 0 ) and the slope intercept of the rain particle size distribution ( ), the prior study correlation shows near independence between these parameters (Figure 2a), whereas the new data indicate significant negative correlation (−0.44) during the convective phase (Figure 2b), when the warm rain processes were dominant, and near independence for the remaining phases (Figures 2c and 2d). This is consistent with the findings presented in Posselt and Bishop (2012), who found that parameter distributions change with time, and points to the need to account for the time-dependent joint parameter distributions in stochastic parameter perturbation schemes. Additionally the parameter variances using the phase-dependent inversions were found to be about one order of magnitude smaller than the variances based on the entire storm evolution (not shown). This provided additional justification for reducing the latter variances to 20% of the original values for the application in this study, as discussed in Section 2.2.
The configurations of the SP scheme considered for the sensitivity experiments are summarized in Table 2. For all configurations the updates of the perturbed parameter values at every time step during the 1D microphysics model simulations were constrained by the boundary conditions for each parameter and the physical relation > . These were imposed by the method of rejection and resampling: invalid samples were rejected and random resampling was applied till the constraints were satisfied. For the configurations using the correlated innovations the resampling involved drawing the vector samples.

Results
The impacts of different configurations of the SP scheme were assessed in terms of the ensemble properties in the observation space and for the state quantities of the 1D column model. The output observable quantities included Accumulated Precipitation (AP), PR, LWP, IWP, OSW and OLW. Note that the AP was not used as an observable in MCMC, but we include it here as it is an outcome of interest for applications; convective storm impact is often as much a function of accumulated precipitation as precipitation rate. The model state quantities included vertical profiles of cloud droplets, rain, pristine ice, snow, and graupel. The impacts in the observation space were measured in terms of time dependent root mean squared error (RMSE) and bias relative to the control simulation and ensemble standard deviation (the spread). The analysis of impacts on the model state quantities was focused on the bias. For all quantities the bias diagnostic was defined as the ensemble mean minus the control (Figure 1).
We note that within the idealized ensemble modeling framework in this study the bias relative to the control simulation represents a mean nonlinear response by the microphysics model, including the transformations into the observations space, to the unbiased stochastic scheme configurations as presented in Section 2.2. This bias measure was used to estimate impacts of the parameter variability represented by the different sampling distributions on the ensemble mean errors in the absence of other factors influencing the errors. The results pertaining to the experiments with 7 perturbed parameters and the decorrelation time scale of 3 hr were selected for the presentation, because the equivalent for the 8-parameter configurations and other decorrelation times were similar. The few notable differences are mentioned in the summary.

Ensemble Performance in Observation Space
The ensemble diagnostics in the observation space exhibited significant time variability for all SP configurations and high sensitivity to the configurations (Figure 3). The time variability was dependent on storm evolution phases and was most pronounced for the PR and radiation observations (Figures 3b, 3e and 3f). For all configurations the bias and spread for the PR exhibited peaks during the early convective phase and then remained high throughout the stratiform phase. The largest PR bias was associated with the SP-LN configuration. For the OLW the bias and spread were small during the early convective phase, large during the transition phase and moderate after (Figure 3e). The highest bias of 8 Wm 2 was produced by the SP-LN configuration. Most storm-phase dependent variability in the ensemble diagnostics was found for the OSW with prominent peaks in the bias and spread during the convective, mid-transition and end-stratiform phases (Figure 3f). The maximum in negative bias during the convective phase was associated with a negative bias in liquid water path, while the peak in negative bias during the transition phase was associated with a negative bias in the IWP (Figure 3c). As for the OLW the highest bias amplitudes were produced by the SP-LN configurations (reaching a minimum of −41 Wm 2 during the early convective phase, Figure 3f). For the IWP, LWP and AP the ensemble spread exhibited steady growth throughout the integration period for all SP configurations, unlike the bias which peaked during the transition phase.
The sensitivity of the ensemble performance to SP configurations was found to be high for all ensemble performance measures and observation quantities. The ensemble bias was overall highest for the SP-LN and SP-N configurations and smallest for the SP-NCovV, while the SP-NCov configuration performance was in between for all observations. The results for the ensemble spread were more dependent on the observation quantity. The SP-NCov and SP-NCovV configurations showed the highest spread values for the AP and IWP, whereas the SP-N and SP-LN configurations were compatible or slightly better for the other observations.
In contrast to the SP configurations the non-SP experiment showed little dependence of the ensemble performance measures on the storm evolution phases (Figure 3, gray curves). As would be expected by the design of this configuration, the ensemble bias was very small and nearly constant throughout the simulation period, with the exception of the last 30 min where it exhibited a slight growth. The RMSE and spread, which were consequently nearly identical, exhibited steady growth throughout the entire period. The spread was significantly smaller than for the SP experiments until about 125 min, after which it reached compatible values and even exceeded the SP-LN and SP-N configurations for the accumulated precipitation. The results point to the fact that the influence of one or more parameters on the model changes with time during the model evolution as various processes related to those parameters are more or less active, leading to a small sensitivity to storm evolution phases when the parameter perturbations are time invariant. The results also suggest that a higher frequency evaluation (e.g., minutes to hourly) of convective scale ensemble prediction is desirable to achieve better understanding of the differences between the microphysics parameter uncertainty representations.
To better understand the relationships between the ensemble properties and the SP configurations the probability density distributions of perturbed parameter values generated in the course of each experiment were examined. The samples consisted of 100 × 2161 realizations based on all ensemble members (100) and all 5-s time steps (2161) for each experiment. For easier comparison between different experiments the marginal probability densities shown in Figure 4 were estimated using the Gaussian KDE (Kernel Density Estimation) with standard deviation value of 4. The primary difference between these exit distributions was in the densities near one or both boundary values. The configurations with the high densities at the extremes of the range such as the SP-LN, SP-N tended to produce a highly biased model response (Figure 3). The exit distributions associated with the SP-NCov configuration exhibited similar properties to the SP-N for most parameters, except and for which the SP-NCov distributions were more symmetrical and with lower densities near the boundaries (Figure 4 blue curves).
In contrast the exit marginal distributions associated with the SP-NCovV configuration which led to the least biased model response had low densities near the range boundaries for all parameters (Figure 4 orange curves). The results suggest that sampling strategies without an imposed covariance structure may lead to more frequent sampling near the extremes, and as a consequence may tend toward a stronger nonlinear biased response. The use of process-level covariance estimates may mitigate this behavior and reduce the incidence of bias.
The results also suggest an inverse relationship for the ensemble spread: high densities near the extremes and mutual independence of parameter perturbations (SP-LN and SP-N configurations) may produce less spread than using the process-dependent co-varied parameters for the same choice of the prescribed entry statistics (the mean and variance values). Stated another way, a restricted range of parameter values produced by the imposition of inter-parameter covariance structures in our experiments led to a larger ensemble spread. We note that the most influential parameters are those that are associated with the snow and graupel fall speeds, and that the coefficients and exponents are strongly correlated (Figure 2). We suspect that ignoring these correlations allows the coefficients and exponents to vary in such a way as to lead to compensating effects, reducing the spread of the ensemble.
The exit distributions also point to the differences between the configurations regarding the impact of the physical constraint < on the sampling. For example, it is clearly evident that using this constraint caused pronounced left and right skewness, respectively, for and distributions for the SP-N configuration (respectively panels a and c in Figure 4, red curves) relative to the Normal distribution that would result from sampling by the linear Gaussian stochastic model according to the expression (Equation 1). As a consequence the parameter perturbations for had a notable bias (shown in Figure 4c by a displacement of red diamond and circle symbols from the control/theoretical-mean marked by black star symbol). Similar but less pronounced tendency for skewed sampling of these parameters was exhibited by the SP-NCov configuration and also SP-NCovV to a lesser degree (Figures 4a and 4c, respectively blue and orange curves and the corresponding mean values). For the SP-LN configuration the impact of the physical constraint was reflected in enhanced skewness for and reduced for . The former was associated with a significant negative bias (Figure 4a, displacement of green diamond and circle from the control). Considering that applying known physical constraints to parameter perturbations is desirable  Table 1. The SP scheme configurations (Table 2) are represented by different colors as shown in the legend. The displayed densities are for the experiments where 7 (solid curves) and 8 (dashed curves) parameters were perturbed using decorrelation time of 3 hr. The theoretical and sample mean values for each parameter are marked by colored symbols along the horizontal axis as follows: theoretical mean (black star), 7 perturbed parameters (filled circles), 8 perturbed parameters (filled diamonds). The symbols are staggered in the vertical to avoid overlapping.
the results point to a potential for mitigating their undesirable impacts on the SP sampling (e.g., introducing or enhancing skewness and bias) by using the process-level covariance estimates with built-in constraints.

Microphysics Response to SP Configurations
In this section the sensitivity of the various cloud hydrometeor species to the SP configurations is analyzed in terms of the differences between the ensemble mean and the control simulation (i.e., the bias). In contrast to the vertically integrated quantities shown in Figure 3, the cloud droplet, rain, pristine ice, snow, and graupel concentrations are more directly related to changes in microphysical parameters. Significant sensitivity to the various parameter perturbation configurations was found for all cloud species outputs. Unlike for the observable (column integral) quantities, the sensitivity for each cloud mixing ratio variable had relatively low time variability, and was associated primarily with the prevalence or not of the specific hydrometeor type in each phase of the simulated storm. For all outputs, with the exception of snow, the biases associated with the different configurations had the same sign and similar spatio-temporal patterns but different amplitudes. Since the microphysical parameters have the most significant and direct effect on the precipitating hydrometeors (rain, snow, and graupel), we analyze these first, followed by the non-precipitating cloud species (cloud droplets and pristine ice crystals).
The bias for the precipitating hydrometeors exhibited the highest variability among the configurations (Figures 5-7). For the rain ( Figure 5) the positive bias (more rain than control) during the transition phase from convective to stratiform (90-120 min) was high for the SP-LN and SP-N configurations (respectively Figures 5a  and 5b), moderate for the SP-NCov (Figure 5c) and negligible for the SP-NCovV and non-SP (respectively Figures 5d and 5e). This positive bias in SP-LN and SP-N configurations is consistent with the high bias in precipitation rate and liquid water path for these two experiments during this time period (respectively Figures 3b  and 3d). The negative rain bias during the stratiform phase had the highest values associated with the SP-LN and lowest for the SP-NCovV, again consistent with the biases (or lack thereof) in PR and LWP in the last hour of the simulation (Figures 3b and 3d). The impact of all SP configurations on the graupel was mass reduction at heights between 6 and 8 km during the transition phase (Figure 6), where the control had the highest load (Figure 1b). The loss of graupel was coupled with the rain increase for the SP-LN, SP-N and SP-NCov configurations (comparing Figures 6a-6c and 5a-5c) but not for the SP-NCovV (Figures 6 and 5d). The result indicates that the different microphysics processes were affected by the parameter perturbations between the different configurations. This was corroborated by the inspection of the differences between impacts of the configurations on the snow  bias ( Figure 7). As for the other hydrometeors the largest bias was associated with the SP-LN configuration (Figure 7a), however the smallest was exhibited by the SP-N configuration which also produced a different spatio-temporal pattern than the other configurations. While the bias for SP-LN, SP-Ncov, SP-NCovV and also non-SP experiments was characterized by a systematic vertical shift in the snow mass profile resulting in the dipole pattern (Figures 7a, 7c-7f), the SP-N experiment produced only a small systematic reduction of the snow mass at approximately the same height during the stratiform phase of the simulated storm (Figure 7b).
The variability of bias amplitudes for the non-precipitating ice and liquid between the SP configurations was smaller than for the precipitating hydrometeors. For the pristine ice the largest impacts were found during the convective phase where the associated microphysical processes were most active (Figure 8). For all configurations the bias was negative for the ice aloft, above about 8 km, and mostly positive below that. This is consistent with the high biases in precipitating hydrometeor contents below; less water vapor mass was available to be transported vertically and converted to pristine ice. Similarly, the impacts of the SP scheme on the cloud droplets for all configurations were significant only during the convective phase (Figure 9). The bias for cloud droplets was uniformly negative and of smaller relative amplitude than for other hydrometeors. Negative bias for cloud droplets was also consistent with the positive bias in precipitating hydrometeors. As for the precipitating hydrometeors the largest bias for the pristine ice and cloud droplets was found for the SP-LN configuration.
As would be expected the impacts of the SP configurations on the ensemble mean microphysics states had a close relationship to the ensemble biases in the observable states. For example, the negative peaks in the OSW bias (Figure 2f) during the convective and transition phases were directly related to the negative biases in the pristine ice and then graupel. Also, the mean loss of graupel mass was reflected in the high negative bias in the IWP (Figure 2d) and the positive in the OLW (Figure 2e). Similarly, the elevation of the snow mass (exhibited as the dipole pattern bias in Figure 7) caused the change of sign for the biases in the IWP and OLW (respectively Figures 2c and 2f) for the SP-LN and SP-NCov configurations which exhibited the significant bias amplitudes for that hydrometeor (Figures 8a and 8c). On the other hand the effect of the non-dipole and low amplitude snow bias for the SP-N configuration during the stratiform phase (Figure 8b) resulted in the low bias for the OSW (Figure 2f, red curves). Overall the results point to the complexity of understanding and ranking the configurations by the mean ensemble properties in the observable quantities alone. Similar to the results in the previous section they recommend using higher frequency evaluations (e.g., minutes to hourly) and also more diverse microphysics-sensitive observations (e.g., space and ground based radar observations) to enable a more accurate capture of different phases of the microphysics in order to better evaluate and understand the differences between the parameter uncertainty representations.

Summary and Conclusions
This study investigates sensitivity of the microphysics parameterization to configurations of the Stochastic Parameters (SP) schemes at the process level. For this purpose we use the generalized linear Gaussian autoregressive order-one stochastic model to represent the stochastic variability in the SP schemes. This formulation allows for exploration of the different configurations of the scheme, including: different parameter distributions, uncorrelated and correlated parameter perturbations, and different decorrelation times within the same algorithm. It is shown that the stochastic process models that are employed in the SPP RP (McCabe et al., 2016) schemes commonly used in the global and regional ensemble systems in prior studies could be represented as special cases of the formulation used in this study. In the current study the stochastic model constitutes the SP scheme.
We apply the scheme to multiple microphysical parameters within the Lagrangian single column version (Posselt & Vukicevic, 2010) of the NASA Goddard Cumulus Ensemble (GCE) model (Tao et al., 2014) for an idealized representation of mid-latitude squall-line convection. This modeling framework was selected because it enables an investigation of the effect of changes in the microphysics parameters on the model output in isolation from any feedback to the cloud-scale dynamics. Also we were able to make use of the prior study results where the same model and test case were used to quantify multi-parameter uncertainty by means of multivariate MCMC inversions using simulated satellite observations of PR, LWP, IWP, OSW and OLW (Posselt & Vukicevic, 2010). This enabled using the estimates of multi-parameter distributions from the prior inversions as the basis for setting the statistical properties of parameter distributions used in the SP scheme. Additionally, the impacts of the non-stochastic and stochastic multi-parameter perturbations on the microphysics model solution and the related satellite-based observable quantities could be directly compared.
The sensitivity to the SP configurations using the 1D column model was carried out through ensemble sensitivity experiments described in Section 2.2 and summarized in Table 2. The experiments explored impacts of using Log-Normal and Normal distributions, uncorrelated and correlated parameter perturbations, number of perturbed parameters, and different decorrelation time scales. The comparison of correlated versus uncorrelated perturbations was performed using only the Normal distribution. The experiments using the correlated perturbations included time-invariant and time-varying covariance estimates. The constant in time covariance was computed using the prior study inversion results, which considered the entire storm evolution. The time-varying covariance was estimated in this study by performing new MCMC inversions separately for the convective, transitional and stratiform phases of the simulated storm. For all experiments the theoretical mean value for each parameter was set to the constant parameter value used in the control version of the deterministic microphysics parameterization, the version that was used to simulate the observations for the inversions with the satellite observations. All configurations were applied to sets of 7 and 8 microphysical parameters and decorrelation times of 1, 3, 6 and 12 hr ( Table 2). For all SP configurations the experiments involved 100-member ensemble simulations where each member was evolved with a different stochastic sequence of parameter perturbations. Additionally the same size ensemble simulations with randomly sampled constant-in-time parameters based on the empirical multi-parameter distribution from the prior study inversion were performed.
The impacts of different configurations of the SP scheme were measured with respect to the control simulation in terms of the ensemble RMSE, bias, and spread in the observation space, and in terms of ensemble mean bias for the the time evolving vertical profiles of rain, graupel, snow, pristine ice and cloud droplets.
The main findings are as follows: • The ensemble diagnostics in the observation space exhibited high storm-phase-dependent time variability for all SP configurations and strong sensitivity to the choice of configuration. • For all SP configurations, the temporal variability was most pronounced for the PR, OSW and OLW observation variables. The peaks in the bias (with respect to the control simulation) and spread for these observables were prominent during different phases of the storm depending on what microphysics processes were most impactful on the particular observable. For example, the bias and spread for OLW exhibited strong increases during the convective phase and then declined through the transition phase. For the IWP, LWP and AP the bias showed pronounced peaks during the convective-to-stratiform transition and stratiform phases, influenced by the changes in the ice microphysics. The ensemble spread for these observables, however, was characterized by low temporal variability and growth throughout the integration period. • Concerning the sensitivity of the ensemble performance to the choice of SP configuration, the most influential factors were found to be the choice of the parameter distribution (normal vs. lognormal) and whether or not the parameter perturbations were correlated. Choice of temporal decorrelation time had little effect on the results of our experiments. For all observables the biases were overall largest for the configurations that used the univariate Log-Normal and Normal distributions (SP-LN and SP-N experiments) and smallest for the configurations with the storm-phase-dependent covariances (the SP-NCovV experiments). For the ensemble spread, the degree of sensitivity to configurations exhibited more dependence on the observable. The SP-NCov and SP-NCovV configurations exhibited the highest spread for the AP and IWP, whereas the SP-N and SP-LN configurations produced compatible or slightly higher spreads for the other observables. • In terms of the RMSE-to-spread performance measure, the SP-NCovV experiments produced the most favorable ensemble response due to least bias and considerable spread. By the same measure the SP-LN configurations provided the least favorable ensembles. • The experiments using constant in time parameters (the non-SP configurations) exhibited very different behavior from the SP experiments: for all observables the ensemble spread had slow growth throughout the simulations and was significantly smaller than for the SP experiments for the first 2 hr. The ensemble bias was small by design because the parameter samples were drawn from a joint empirical parameter distribution constrained by the simulated observations over the entire storm evolution and centered on the parameter mean values. The results suggest that time-invariant parameter perturbations do not support the variability of model responses due to different microphysics processes during the storm evolution. • Significant sensitivity to the configurations was found for all microphysics state outputs. Unlike for the column integral observable quantities, the sensitivity for each cloud microphysical mass mixing ratio variable had relatively low time variability, and was associated primarily with the prevalence or not of the specific hydrometeor type in each phase of the simulated storm. For all outputs with the exception of the snow, the biases (with respect to the control simulation) associated with the different configurations had the same sign and similar spatio-temporal patterns but different amplitudes. • Overall, the bias in the microphysical states was highest for the SP-LN and SP-N configurations and lowest for SP-NCov and SP-NCovV. Most significant biases were exhibited for the rain, graupel and snow. • The sensitivity to changing the number of parameters from 7 to 8 by adding the slope intercept for the rain particle size distribution was negligible overall except it caused higher time variability of the ensemble spread and RMSE for the precipitation rate.
These findings lead to the following main conclusions and recommendations.
• The impact of the stochastic parameter perturbations on the microphysics model ensemble solution is highly storm-phase dependent • The correlated parameter perturbations informed by the process-level covariance estimates may reduce biased response by the microphysics parameterization and introduce more variability than the uncorrelated sampling using the univariate distributions • Using the Log-Normal distribution or other skewed sampling may lead to undesirable biased response by the microphysics parameterization • Using the appropriate mean and covariance estimates for the stochastic parameter sampling for the different phases of the microphysics dependent on the storm morphology is highly desirable • This calls for the parameter estimation capabilities using the microphysics-sensitive observations preferably with hourly or less temporal frequency such as geostationary satellite and ground-based radar remote sensing. Such capability could be enabled by convective scale data assimilation systems The primary caveat to this study, and the factor limiting its applicability in real-world prediction scenarios, is the use of a single column model framework. While this has enabled us to conduct a wide range of experiments and has enabled a more straightforward interpretation of the results, there is no feedback from parameter perturbations to the three-dimensional dynamics and thermodynamic state. This is also likely the reason there was little to no sensitivity of our results to changes in the decorrelation time scale; the 1D model evolution is controlled by the imposed water vapor source and vertical velocity profile. In previous parameter estimation studies, we have shown that the parameters that are most influential, the time during which they have an effect, and the inter-parameter covariances are quite similar between 1D and 3D models (compare results reported by Posselt and Vukicevic (2010) with Posselt (2016)). However, these studies considered only non-stochastic perturbation methods. The natural next step for our study is to employ the various SP methodologies in a three dimensional simulation framework and perhaps also to examine the ensemble spread due to initial conditions versus parameter perturbations as was done for static parameter perturbations by Posselt et al. (2019). We expect that the temporal autocorrelation length scale will have a much more substantial effect in a model that allows feedback from the parameterization to the thermodynamic and dynamic environment.
In addition, we have examined microphysical parameter perturbations in one particular type of weather system (deep convection). It would be useful to explore the influence of choices in SP scheme for other weather systems (e.g., extratropical cyclones and/or orographic precipitation). We also are aware that cloud microphysical schemes are not the only parameterizations with uncertainty. Other parts of the model (e.g., the planetary boundary layer parameterization; Aksoy et al., 2006) contain uncertain parameters and also would benefit from study of various SP scheme configurations. We plan to pursue these and other studies in the near future.

Appendix A: Relationship Between Generalized Linear Gaussian AR1 Formulation and the SPP and PR Schemes
In this appendix the correspondence between the SPP and RP schemes and the generalized linear Gaussian AR1 formulation used in this study is presented. By the generalized formulation the stochastic parameter perturbations are evolved according to the expression (Equation 1), repeated here for convenience where Φ is a vector stochastic quantity with mean Φ , is the first order regression coefficient determined by the model time step and decorrelation time scale Λ , 1∕2 Φ ( ) is the innovation term where 1∕2 Φ is square root covariance matrix, and ( ) is a (0, I) random vector, with being the identity matrix. Assuming no spatial variability for the parameter perturbations the SPP scheme formulation would derive directly from the expression (Equation A1) by using a diagonal variance matrix Φ and = − Λ . As discussed in Section 2.2, the definition of the stochastic quantity would depend on the choice of statistical distribution for the parameter variability within a range of permissible values. For the scheme application in the ECMWF medium-range ensemble forecasting system Ollinaho et al. (2017) used the Long Normal distribution. For this distribution as is shown in Section 2.2 the vector Φ consists of the quantities = ( ′ ) , where and ′ are, respectively, the unperturbed and perturbed parameter values for each parameter to be perturbed. Ollinaho et al. (2017) considered two options for based on the following two criteria: ′ = and ′ = . The first criterion corresponds to the SP-LN configuration in this study with = − 2 2 and 2 = ( 2 2 + 1 ) , where 2 is predefined parameter variance, whereas the second criterion corresponds to = 0 .
The recent applications of the SPP scheme with the regional models (Jankov et al., 2019;Stanford et al., 2019;Thompson et al., 2021;Wang et al., 2019) considered the Normal distribution for the SPP using the multiplicative stochastic noise formulation for the relationship between the perturbed and unperturbed parameter values according to ′ = (1 + ) . For this formulation the expression (Equation A1) would apply to = ′ − 1 , = 0 and = , equivalent to the SP-N configuration used in this study. When the spatial variability of parameter perturbations is considered and represented by the Gaussian correlated noise pattern in 2D horizontal space of the prediction model, as was done for all prior applications of the SPP scheme, the stochastic quantity would involve spectral coefficients of the spatial noise pattern.
The PR scheme (Hermoso et al., 2021;McCabe et al., 2016) assumes no spatial variability of the parameters and makes use of the uniform distribution for the stochastic variability according to ( + 1) = ( ) + 2 √ 1 − 2 where = 1 − Λ and is a random quantity with uniform distribution in the interval [−1,1], could be represented by the model (Equation A1) using = ′ − 1 with mean = 0 . To represent the sampling equivalent to using the bounded uniform distribution the variance may be set to = 2 * √ 1 − 2 , where * would be a large value applicable to all parameters. The RP scheme and this configuration are not exact equivalents due to approximating the uniform distribution by the Normal with a large variance, but the correspondence would be close for all practical purposes assuming application of the appropriate boundary conditions (the range of parameter values). The RP configuration was not explicitly considered in the sensitivity experiments in this study because it is similar to the SP-N configuration (Table 2).

Data Availability Statement
The data used in creation of this manuscript include: simulations by the 1D column microphysics model used in this study, diagnostics derived from the model simulations and results from Bayesian inversions performed in Posselt and Vukicevic (2010) and in this study using the same model and simulated column integral satellite observations. The data are available at https://doi.org/10.5281/zenodo.5736940. A guide to content and format of the archived data is provided in the Supporting information. The graphics were generated using Matplotlib (https://matplotlib.org) and Seaborn (https://seaborn.pydata.org) software libraries. The 1D column microphysics model used in this study is a single column version of the NASA Goddard Cumulus Ensemble (GCE) model (Tao et al., 2014), described in Posselt and Vukicevic (2010). The Bayesian inversion model used in this study is presented in Posselt and Vukicevic (2010).