Numerically Relevant Timescales in the MG2 Microphysics Model

Climate models rely on parameterizations of a variety of processes in the atmospheric physics, but a common concern is that the temporal resolution is too coarse to consistently resolve the behavior that individual parameterizations are designed to capture. This study examines timescales numerically derived from the Morrison‐Gettelman (MG2) microphysics as implemented within the Energy Exascale Earth System Model, version 1 (E3SMv1). Numerically relevant timescales in MG2 are derived by computing the eigenspectrum of its Jacobian. These timescales are found to often be smaller than the default 5 min time step used for MG2. The fast timescales are then heuristically connected to individual microphysics processes. By substepping a few particular rain processes within MG2, the time discretization error for those processes was considerably reduced with minimal additional expense to the overall microphysics. While this improvement has a substantial effect on the target processes and on the vertical distribution of stratiform‐derived rain within E3SMv1, the overall model climate is found to not be sensitive to MG2 time step. We hypothesize that this is because the surface climate does not depend strongly on certain process rates, especially MG2's rain evaporation rate.


Introduction
Atmospheric general circulation models (GCMs) consist of resolved-scale fluid "dynamics" and parameterized "physics" components. Dynamics is based on the Navier-Stokes equations, which makes its discretization and numerical solution straightforward (though computationally taxing). In particular, the dynamics time steps are typically controlled by the Courant-Friedrichs-Lewy condition (Dennis et al., 2012) and other conditions that are needed for stability at a given vertical and horizontal resolution (e.g., the need within a semi-Lagrangian advection scheme to avoid evacuating all mass from a layer). Parameterized physics, on the other hand, handles the grid-scale effects of the remaining collection of subgrid-scale atmospheric processes. Because the proper equations for these processes are often not known, and because the processes themselves often do not behave smoothly (e.g., across the clear-sky/condensate boundary), numerical treatment of time integration in atmospheric physics has received much less attention. Much of the physics time step analysis that has been done focuses on issues of dynamics-physics coupling. This may be either in an idealized setting (Dubal et al., 2004;2005;2006;Staniforth et al., 2002a;2002b) or using full-complexity models, where the coupling between model dynamics and convection parameterizations in particular can cause models to become sensitive to the time step (Mishra et al., 2008;Mishra & Shanay, 2011; 10.1029/2019MS001972 Williamson & Olson, 2003;Williamson, 2013). That said, some work has also been done on process coupling within the physics specifically (Beljaars et al., 2004;Wan et al., 2015). This work often focuses on the effects of using a sequential split approach, where the model state is updated by each physics parameterization in turn before the next parameterization acts; for instance, when using this form of coupling, the order in which parameterizations are evaluated can have a significant effect on results (Beljaars, 1991;Donahue & Caldwell, 2018).
The time steps used for individual parameterizations within a model's physics are less likely to be chosen based on stability considerations compared to the dynamics. One reason is that limiters can prevent theoretically unstable parameterizations from producing unrealistic average outputs even at arbitrarily long time steps (though some numerical artifacts may remain). For instance, Morrison and Gettelman (2008) showed that using a long time step for a cloud microphysics scheme results in instability that causes wild oscillations in the precipitation rate over time, but with almost no change in the mean. As a result, physics time steps are in practice often set based on a tradeoff between throughput requirements and sensitivity of the mean values to further time step reductions, rather than by focusing on the stability of a given parameterization.
A more formal numerical analysis of atmospheric physics routines would be of great interest. Traditional stability analysis can provide guidance regarding the time steps that should be used for various parameterizations. When solving a linear problem, we can readily determine whether any given method is absolutely stable (i.e., when the error does not grow exponentially in time), by checking the eigenvalues associated with this problem. In the nonlinear case, there is no simple procedure that rigorously bounds the error growth in this way, but if the equations governing process rates are smooth, it can still be informative to analyze stability using the linearization of the problem. In this study, we will make heavy use of the fact that the reciprocals of the eigenvalues of the Jacobian matrix are key timescales for any linearized dynamical system, and are typically good estimates for the timescales of the original nonlinear problem as well (LeVeque, 2007). Using a time step size small enough to solve the linearized system will not guarantee any particular level of accuracy, but we can regard it as a necessary condition. If the numerical method being used would ordinarily be unstable at a given time step, and model stability is only being maintained by limiters to prevent this instability, it is likely that the model's time step is too large to accurately approximate the solution of the underlying mathematical equations for the parameterized physics.
The current study uses this approach to analyze version 2 of the Morrison-Gettelman microphysics (MG2) , the component of the Energy Exascale Earth System Model version 1 (E3SMv1) responsible for microphysics in stratiform clouds. Previous work on the Community Atmosphere Model version 5 (CAM5), a predecessor to E3SMv1, shows that cloud cover and cloud ice distribution are significantly affected by the model time step (Wan et al., 2014) and also that the stratiform microphysics considerably affects the temporal convergence rate of the overall model (Wan et al., 2015). CAM5 used an earlier version of the Morrison-Gettelman microphysics (MG1), which used a 15 min time step size. This step size was originally chosen because important quantities (e.g., radiative cloud forcing, liquid water path) appeared to have converged in the mean . However, this was based on a fairly modest range of time step sizes (either one, two, or three microphysics substeps per 30 min model time step).
When MG2 was introduced to replace the original MG1 microphysics, it added prognostic precipitation, and the model time step was reduced to 300 s in both E3SMv1 and CAM5s successor model, CAM6. In , tests of MG2 in an idealized test driver suggested a sensitivity to the vertical distribution of cloud and rain mass to time step size when the step size is 120 s or greater, as well as a significant dependence of liquid water path on time step, though the authors attribute this to the simplified condensation present in their driver. In , the authors note some benefits of reducing the time step of both the microphysics and macrophysics parameterizations to 300 s when using MG2 in CAM5, at the cost of a 21% increase in computational cost for that model. This trade-off between accuracy and affordability motivated the use of a shorter microphysics time step.
There are relatively few studies on the time step demands of bulk microphysics schemes with prognostic precipitation, especially when time step sizes of several minutes or more are used. Older models that ran at longer time steps tended to use diagnostic precipitation as MG1 did, though there are exceptions (Lopez, 2002). In most cases when prognostic precipitation is used, precipitation-related processes are substepped at a smaller time step, no more than 120 s (Fowler et al., 1996;Michibata et al., 2019;Sant et al., 2015).

10.1029/2019MS001972
In Posselt and Lohmann (2008), prognostic rain was added to the ECHAM5 model, and this was tested at time steps up to 300 s, finding that a time step of 90 s or less was necessary to reach the converged rates of autoconversion and accretion. Another interesting result can be found in Chosson et al. (2014), which found a 60 s time step to be an adequate time step for a scheme with prognostic precipitation, but found a 120 s time step to be inadequate.
Note also that MG2 relies on limiters to avoid negative cloud liquid mass. In particular, , which first tested the scheme in CAM5, found that for a mixed-phase single column case, the limiter still triggered at time steps as low as 75 s, due to excessive liquid evaporation and vapor deposition onto ice (as well as some direct accretion of liquid by snow). This also motivates a need for a more detailed investigation of the time step sensitivity of E3SMv1's cloud processes, since the limiter should rarely have a significant effect if the microphysics is adequately resolved.
In this paper, section 2 will describe the relevant features of MG2 and its usage within E3SMv1. Then we will address questions about the numerically relevant timescales in MG2 in several sections. First, we examine the eigenvalues of the Jacobian in section 3, which are derived numerically based on a broad sample of conditions from the global model. Then, we associate these eigenvalues with specific processes in section 4. This allows us to connect the timescales associated with MG2 to subsets of the specific physical processes it implements. Next, in section 5 we examine the timescales of MG2 in different regimes where specific processes are dominant (e.g., warm versus cool grid cells, cloudy grid cells versus precipitation in an otherwise clear sky). Finally, in section 6, we check these results by substepping the processes identified as having fast timescales. Each of these four sections is divided between one or more methodology sections, followed by a section presenting the results. We will conclude by discussing the relevance of these results to the numerics of MG2 generally and introduce some preliminary information about the impact of MG2's time step on the global model.

Model Description
E3SM is an ongoing U.S. Department of Energy project designed to produce a state-of-the-art Earth system model that can leverage the Department of Energy's largest supercomputers to produce high-resolution simulations. The scientific goals of this project relate to three main topics: (1) the water cycle, (2) the cryosphere, and (3) biogeochemistry (Golaz et al., 2019). While most of this paper is concerned with its stratiform microphysics (MG2) running in isolation using a stand-alone driver, we do use the global model to generate input data and do some preliminary investigation of the impact of MG2's time step on global climate. For this purpose, we use version one of the E3SM Atmosphere Model (EAMv1), using prescribed data for the sea surface temperature and sea ice extent, and running at an approximately 100 km (1 • ) horizontal resolution (Rasch et al., 2019;Xie et al., 2018).
Within EAMv1 at this resolution, much of the physics, including the deep convection scheme, runs at a 30 min time step, which is also the frequency of physics-dynamics coupling. However, both MG2 and the Cloud Layers Unified By Binormals parameterization (known as CLUBB, which handles cloud macrophysics, shallow convection, and turbulence) use a shorter time step. These parameterizations, as well as a handful of aerosol-and ice-related microphysical processes calculated outside of MG2, are substepped in a loop within the physics driver, and use a 5 min (300 s) time step. Although the time step for much of the physics is reduced to 15 min when running EAMv1 at high resolution, the MG2 and CLUBB time step is not changed. We therefore regard this as the default MG2 time step and will concern ourselves with the behavior of MG2 only at 5 min and smaller time step sizes.
The MG2 microphysics uses a two-moment representation of four hydrometeor types. For cloud liquid, the two prognostic variables are the grid cell average mass mixing ratio (q c ) and number concentration (n c ), and process rates are calculated assuming that the particle diameter follows a gamma distribution, with the mean size dictated by the ratio q c ∕n c , and the other shape parameter diagnosed from the number concentration. Similarly, an average mass mixing ratio and number concentration are used to describe the model's cloud ice (q i , n i ), rain (q r , n r ), and snow (q s , n s ), though these other hydrometeors are assumed to have sizes that follow a simpler exponential distribution. The temperature (T), the humidity (q), and these hydrometeor variables together make up the 10 prognostic state variables for MG2, and for the purposes of this paper we define the "state" used by MG2 to consist of only these variables. Other inputs are also used by MG2, such as Size limiters n c , n i , n r , n s Limiters constraining hydrometeor particle sizes to remain in relevant ranges the cloud fraction, but these are diagnosed by other parameterizations, and we assume that they are roughly constant over the course of the 300 s time step at which MG2 runs.
Parallel splitting is used for most of the processes in MG2; that is, the inputs to these processes are generally the same as the inputs to MG2 itself. Once the process rates are calculated, the state is updated by adding contributions from all processes and taking a single forward Euler method step. As part of this process, a series of conservation checks are performed. For each hydrometeor, MG2 checks whether the forward Euler step would produce a negative mass mixing ratio or number concentration (e.g., using up all cloud liquid). If so, it calculates the ratio between the amount of mass (or number) actually present, and the amount that the processes attempted to subtract. The process rates are scaled down by this same ratio, which allows MG2 to avoid producing a negative concentration without changing the relative strength of different processes. All processes that are applied at this first, parallel stage are listed in Table 1. Note that the droplet activation is a special case; its rate is calculated outside of MG2, and it is applied sequentially using a single forward Euler step before all the other processes in this table. There are also three other processes controlled by "external" schemes. While these are applied as if they had been calculated within MG2 itself, within E3SMv1 they are prescribed by a different scheme.
After this, MG2's sedimentation is run on the updated state. The sedimentation calculates the mean fall speed of each hydrometeor in each grid cell and uses a first-order explicit upwind scheme to update the profile of the entire column. Hydrometeors that reach the surface are considered precipitation and reported as a surface flux to the host model.
Some limiters and instantaneous adjustments (e.g., forcing all precipitation to freeze/melt when introduced to very cold/warm grid cells) are also applied at two stages: at the beginning of MG2, and after the sedimentation.

Timescales in MG2
In this section we investigate the inherent timescales present in the MG2 code by calculating the eigenvalues of a numerically derived Jacobian.

Methodology 3.1.1. Running MG2 as a Stand-Alone Process
Within E3SM, MG2 and CLUBB are substepped together using a 300 s time step, which is both the default and maximum recommended time step for CLUBB. We can view the main role of MG2 as providing the average rates of change of MG2's 10 state variables due to microphysical processes over a given time step. We can evaluate how well MG2 answers this question without use of the full E3SM model, as long as we have a representative sample of atmospheric columns and a method for running MG2 outside of E3SM.
To obtain a sample of input columns, we ran a standard preindustrial global E3SM simulation with ∼100 km atmospheric resolution (ne30_ne30 grid) and prescribed sea surface temperature (compset F1850C5AV1C-04P2). The code used for this integration was a beta version (git hash 7a17edbe) which is structurally similar to the E3SMv1 release. Code modifications between this version and the release are not expected to make a difference for the results here. At the end of a 5 day run, we saved a "snapshot" of every column on the planet at the final time step, including every input used in the final call to MG2.
This provided us with a representative sampling of columns for a particular day in January, which is the basis for the analysis in the remainder of this paper. A larger or more varied sample of columns, taken over the course of a run lasting a year or more, would have some statistical differences, especially since such a sample would capture seasonal variation. However, since this sample causes MG2's most important processes to operate under a wide variety of conditions, we do not believe that broader sampling would affect our overall conclusions.
In order to run MG2 in isolation from the rest of the model, we used F2PY to create a Python interface and compiled MG2 and the interface into a stand-alone library. This library was used by a set of drivers to produce the MG2 stand-alone results in this paper.
We narrowed down the number of active processes to a subset of those in MG2 by implementing switches to disable two types of process. First, sedimentation was disabled in our runs. MG2's sedimentation runs sequentially after the state is updated by all other processes, and it uses an adaptive time step, smaller than that of MG2 as a whole, to satisfy its Courant-Friedrichs-Lewy condition. The numerical problems posed by the sedimentation are of a different character from those posed by the rest of MG2, and so we believe that it would be best to examine the sedimentation in a separate study. In fact, we consider the sedimentation to not be a microphysical process at all, but rather a dynamical process which happens to be packaged in MG2.
We do not believe that enabling sedimentation would affect the overall conclusions that we draw about the timescales relevant to resolved microphysics, especially since (as we shall see), these timescales are often too short to allow much sedimentation to occur. However, since the sedimentation time step can still be quite short for heavy rain (<10 s), we do note that there may be benefits to coupling the sedimentation to MG2's other processes more frequently, though evaluating such a change is beyond the scope of this paper.
Since sedimentation is the (only) process in the MG2 microphysics that transports mass between levels, disabling sedimentation allowed us to view MG2 as a tool for solving a set of ordinary differential equations in a collection of uncoupled grid cells, rather than as a tool for solving a set of partial differential equations on a 1-D atmospheric column. This provides two benefits: 1. All remaining processes have a well-defined and easily controlled time step, since only the sedimentation uses an adaptive time step. 2. If sedimentation was enabled with 72 vertical layers per column, then the state vector input to MG2 would be 720-dimensional (given the 10-dimensional state for each grid cell). With sedimentation disabled, MG2 operated on individual grid cells with a 10-dimensional state in each one.
Second, many "instantaneous" processes were disabled. Instantaneous processes are those processes which MG2 implements not by calculating an expected process rate but by making an immediate adjustment from a state that is considered unstable or a violation of model assumptions, to a more stable or valid state. Usually such processes involve rapid changes of state, for example, quickly melting snow that has fallen into a very warm grid cell.
We disabled these processes because we were interested in how the microphysical processes in MG2 are resolved, but instantaneous processes are by definition so fast that the model is not designed to resolve them in the first place. Furthermore, disabling the instantaneous processes removed any explicit dependence of the calculation on the time step size, which allowed the time step chosen during the Jacobian calculation to be chosen arbitrarily without impacting the results. Since these processes generally involve rapid melting/freezing of hydrometeors, we believe that our conclusions for warm cases, without significant ice physics, should not be affected by this change. However, when ice physics is present, MG2's resolved processes may be exposed to some states that are very different from those normally present in an E3SM run (e.g., large amounts of rain in very cold grid cells), so we should note that this modification results in MG2 effectively operating in a larger range of regimes than those normally present in a full model run. We take this behavior into account in interpreting results, so our conclusions are not erroneously affected by this modification.
One final change that was made to MG2 was to extract the calculation of precipitation area fraction from the main code, so that it was only calculated once per driver call. In conjunction with disabling sedimentation, this change was needed to allow MG2 to be run on individual grid cells rather than full columns, which greatly reduces computational cost. Furthermore, MG2's default precipitation area fraction calculation depends mostly on the cloud fraction calculated by CLUBB and is insensitive to changes in the state vector. As a result, forcing this variable to remain constant for short MG2 runs does not produce noticeably different results from allowing it to vary over time.

Measuring Differences in MG2 States
We often wanted to compare MG2 states in order to (a) establish whether MG2 is active enough to produce a final state that is very different from its input state, or (b) measure the magnitude of a change in MG2's output given a change in MG2's parameters. However, MG2's 10 state variables include temperature, mass mixing ratios, and number concentrations, which have different units and orders of magnitude, so an isotropic norm such as the L 2 norm is ill-suited to represent the distance between MG2 states. One way to remedy this would be to use a set of constant weights to convert all variables to a common set of units; for instance, a difference in temperature could be converted to an approximate change in mass mixing ratio of liquid water that would be necessary to produce that change in temperature via evaporation/condensation.
A simpler approach is to focus solely on the mass mixing ratio terms in the state vector, from which we define a quantity called the total water mass difference (D w ). This is simply half of the L 1 -norm of the water mass, or for two states labeled as s 1 and s 2 : (1) This value can be interpreted as the total amount of water that is in a different category between two states. For instance, consider if two states differ due to the evaporation of 10 −3 g/kg of rain. Then q r is reduced by 10 −3 g/kg, while q is increased by the same amount; after dividing by 2, we get D w = 10 −3 g/kg.
The total water mass difference is arguably the most straightforward method for measuring the differences between MG2 outputs. Technically, it does not account directly for differences in temperature, or at all for differences in number concentration. However, MG2 can only create changes in temperature via phase change, which will affect the mass mixing ratios, and in practice changes in number concentration will only occur in circumstances where processes that affect mass mixing ratio are also present. Differences in the overall MG2 state can therefore be distinguished by looking at mass mixing ratio alone.

Jacobian Eigenvalue Calculation
A standard method for analyzing the stability of a numerical method on a nonlinear problem is to linearize the system about a given state by calculating the Jacobian (LeVeque, 2007). For a Jacobian with real eigenvalues, the eigenvectors represent those directions in which a perturbation to a system causes purely positive or negative feedback in the same direction, with the associated eigenvalues representing the sign and strength of that feedback. When an inappropriate numerical method is used, or the time step used with a given method is too long, the system's response to a perturbation will tend to rapidly amplify time integration errors, leading to instability. Accordingly, the stability of a linearized problem depends on the location of the eigenvalues of the Jacobian (multiplied by the time step) in the complex plane, and specifically on whether they fall in the stability region for a particular method (in our case, forward Euler).
Alternatively, one can think of the inverse of these eigenvalues as a set of timescales associated with the linearized system. In order to keep the system stable without depending on limiters or other artificial corrections to a given method, our time step must be sufficiently small compared to the smallest timescales derived from the Jacobian, so that the product of time step and eigenvalue lies within the method's stability region.
However, calculating this Jacobian presents some difficulty. The system of equations represented by MG2 is quite complex, and contains a large number of thresholds that adversely affect the smoothness of the process rates (e.g., the process rates are generally continuous, but not continuously differentiable, at the freezing point of water).
Rather than attempt to analytically differentiate MG2 around each point, we numerically calculated the Jacobian using the numdifftools Python package. The specific method we used from this package is a second-order forward difference method, supplemented with one stage of Richardson extrapolation (so the resulting method is third-order accurate).
The use of a forward (rather than central) difference method was important. MG2 produces floating point exceptions when given negative concentrations, so we needed to use a one-sided method in order to numerically linearize about a state where any constituent has zero mass. Furthermore, by using a small linear change of variables, we could also ensure that this one-sided method does not use perturbed states that violate MG2's size limiters. Since MG2 enforces these limiters at several points throughout the code, any perturbation that leads to a state outside of the acceptable size ranges will result in an instantaneous adjustment back to a valid state, and these instantaneous adjustments are not of interest.
Since the Jacobian is a 10 × 10 matrix, finding its eigenvalues and eigenvectors was a negligible cost compared to the MG2 calls necessary to calculate the Jacobian in the first place. Thus, we simply used SciPy's linalg.eig (which wraps LAPACK calls that use the QR algorithm) for this calculation. The MG2 Jacobian almost always has a full set of eigenvectors when MG2 is active, though the Jacobian becomes defective when no hydrometeors are present. This appears to be because the process rates are completely insensitive to the temperature and humidity variables (rates are always zero in the absence of hydrometeors), but perturbations to the hydrometeor masses will result in temperature and humidity changes. This introduces asymmetric off-diagonal elements to the Jacobian that make it nondiagonalizable. We had no interest, however, in states where MG2 is inactive or nearly inactive. We therefore only included cases where the process rates are sufficiently large, as defined by a total water mass difference threshold. Specifically, for an initial state s 1 and a final state s 2 we were only interested in columns where D w (s 1 , s 2 )∕Δt ≥ 10 −7 g kg −1 s −1 . Since Green bars represent eigenvalues associated with at least one active process, but no "primary" process. Blue bars represent eigenvalues associated with inactive processes. A black line is placed at (300 s) −1 for comparison with the MG2 time step.
most grid cells in the model have little to no condensate, this eliminated most of the total grid cells from consideration (79.7%), which also significantly sped up calculations on the full sample.
The eigenvalues of the Jacobian are furthermore almost always real or have small imaginary parts, as shown in Figure 1. This figure also shows that the few eigenvalues that have a large imaginary part also almost always fall within the region of absolute stability for the forward Euler method at a 300 s time step, unlike many of the purely real eigenvalues, which often have much larger magnitudes. We are therefore solely concerned with the real parts of these eigenvalues from this point forward. Figures 2 and 3 contain histograms of MG2's eigenvalues for grid cells above the 10 −7 g/kg/s cutoff. These eigenvalues are further categorized based on whether they are associated with active or inactive processes, where an active process is one that either affects water mass above a rate of 10 −7 g/kg/s, or affects a hydrometeor number above a rate of 2.98×10 −3 /kg/s, which is equivalent to the rate of 400 μm diameter rain particles Figure 3. Histogram of the real parts of the eigenvalues of MG2's Jacobian, focusing on positive eigenvalues. Red bars represent eigenvalues associated primarily with an active process (see section 4 for further details on this association). Green bars represent eigenvalues associated with at least one active process, but no "primary" process. Blue bars represent eigenvalues associated with inactive processes.

Results
that would have to be introduced to produce a mass change of 10 −7 g/kg/s. The details of how eigenvalues are associated with processes are explained in section 4; we use association here just to indicate which eigenvalues are physically meaningful and which are akin to numerical noise. We can roughly divide these eigenvalues into five categories: 1. Eigenvalues associated with inactive processes (shown in blue). These eigenvalues are typically related to physics that is not active in a given regime, for example, ice physics in a warm grid cell, so they are not as relevant to the numerics in practice. 2. Negative eigenvalues of magnitude greater than (300 s) −1 . These eigenvalues correspond to short-timescale processes, which have rates that may decay too rapidly for the default MG2 timescale to handle. 3. Negative eigenvalues of magnitude less than (300 s) −1 . These eigenvalues correspond to processes which MG2 is able to resolve, assuming roughly linear behavior of MG2. 4. Near-zero eigenvalues. These eigenvalues correspond either to extremely slow feedbacks within MG2, or to forbidden directions of motion in the phase space (particularly eigenvectors that are perpendicular to surfaces of constant energy or mass, since MG2's process rates must be tangent to such surfaces). 5. Positive eigenvalues. These eigenvalues correspond to eigenvalues of processes that are temporarily in a state of positive feedback (e.g., accretion produces larger raindrops, and those larger raindrops will be effective at accreting even more cloud water). Generally speaking, MG2 avoids instability from these processes due to its nonlinearity, since all MG2 processes "use up" some form of mass or number, and therefore will eventually slow down over time.
In general we were most interested in dealing with the second case, the large negative eigenvalues, since these are common and correspond to cases where MG2's large timescale is likely to cause problems. Note that while the forward Euler method is absolutely stable up to eigenvalues of (150 s) −1 at this time step, the (300 s) −1 threshold is the point at which the method will produce a monotonic decay in the error. This is more relevant for avoiding negative mass concentrations, and therefore for avoiding triggering limiters that artificially constrain the process rates. Figure 2 shows that there are a great number of timescales that are not adequately resolved at MG's current time step (by this measure, even more than the number that do seem adequately resolved).
We also wanted to understand the nature of the positive eigenvalues, since these eigenvalues are always outside the region of stability for the forward Euler method (though this may not be as bad as it seems, as we discuss in section 4.2).
In order to better interpret this result, however, we needed to associate the eigenvalues more closely with both specific sets of processes within MG2, and the physical conditions under which these eigenvalues arise. For instance, there may be cases where the model is formally unstable without limiters at a given time step, but in practice there is little difference between the resolved and limited behaviors.

Connecting Timescales to Specific Processes
Documenting the timescales of MG2 processes was of inherent interest, but we were particularly interested in the finding that MG2 is often integrated using a Δt which is too large to represent the processes it represents. Substepping all of MG2 in order to capture these timescales is not computationally feasible, but numerically accurate solutions may still be affordable if the need for substepping could be isolated to just a few processes.

Methodology 4.1.1. Measuring Eigenvalue-to-Process Associations
We sought to assign each eigenvalue from section 3 to one or more of the MG2 processes listed in Table 1. Let us label the grid cell output tendencies from applying MG2 to state s as r = s∕ t. The Jacobian J r has i, jth entry r i ∕ s j where r i is the ith entry in r and s j is the jth entry in s. As noted above, J r is diagonalizable for the states that were of interest for us, so matrices of eigenvalues Λ and eigenvectors V exist such that VΛV −1 = J r .
If we label the tendencies due to a particular process p as r p , then Σ P p=1 r p = r. This leads directly to the identity 10.1029/2019MS001972 which is the heart of our association method. Now construct a matrixC of dimensions N (number of eigenvalues) by P (number of processes to consider) whose n, pth entryC np is the nth diagonal element of V −1 (J r p )V. Since it depends only on these diagonal elements,C is independent of the scaling of the columns of V. Additionally, so the nth row ofC represents the contributions from each process to the value of n . The elements that contribute to a particular eigenvalue may be of different signs. Two large-magnitude values ofC np may add to produce a larger n , or they may partially cancel to produce a smaller n . In either case, however, we would say that the largest magnitude values inC have the most influence on the corresponding eigenvalues.
As a final step, normalizeC by taking the absolute value of each element, and by dividing each row by its 1-norm, producing a new matrix C. Then each element of C describes the fractional contribution of its column's process to its row's eigenvalue. If any element of C is greater than 0.5 (our heuristic for whether at least 50% of an eigenvalue's magnitude can be attributed to a particular process), then we say that the eigenvalue for that row is primarily associated with the process for that column. This is how the colors in Figures 2 and 3 are derived: the red bars count eigenvalues primarily associated with active processes, blue bars correspond to inactive processes, and green bars represent eigenvalues without a primary association.
We are mainly concerned with large-magnitude eigenvalues, which almost all have a primary association with a particular process.
For purposes of this study, we treated the effects of MG2's size limiters as if they were a separate physical process. This was largely because MG2's state can become unrealistic or invalid if these limiters are disabled, so we were not able to disable them as we had with other instantaneous processes, and it was therefore necessary to account for their effects on MG2's state. Note also that these limiters are not applied purely for artificial, numerical reasons, since the maximum size limiters are also the means by which MG2 accounts for the spontaneous breakup of large particles, which is important especially for a reasonable treatment of precipitation.

Correlation of Multiple Processes
In addition to identifying processes associated with particular timescales, we wanted to identify tightly coupled processes. Such processes must be solved simultaneously or substepped together in order to obtain numerically accurate solutions. In addition, accounting for such coupling is needed in order to create useful conceptual models of microphysical behavior.
We hypothesized that processes that are tightly coupled in this way would be tend to be associated with the same eigenvalues. That is, while most eigenvalues have a single primary association with another process, there are also large values of C linking these eigenvalues to other processes, and we believed that this might signify cases where both processes must be solved together for accurate time integration.
To identify such possible sets of tightly coupled processes, we therefore first tallied the number of eigenvalues with primary associations to each process. We then looked at the average value of C across all eigenvalues primarily associated with a given process. If, for instance, eigenvalues that are primarily associated with autoconversion also typically have strong association with accretion, then we would focus on autoconversion/accretion coupling for further study and optimization.

Results
We will first examine the association between positive eigenvalues and specific processes, shown in Figure 4.
First, we found positive eigenvalues that are associated mainly with accretion-related processes (e.g., note the large number of eigenvalues associated with liquid accumulation onto rain and snow in Figure 4). These eigenvalues are positive when the accumulating particles are small and few in number, since in this case the initial accretion increases the size of the particles, making them better at accumulating additional mass.
In the long run, this is not a threat to the stability of the model, since eventually the accretion will begin to deplete the cloud, which will cause the rate of accretion to slow (see Figure 5). We will not examine this issue further here, but it should be noted that long time steps may delay the onset of heavy precipitation in these cases. Second, we found eigenvalues that are associated mainly with evaporation or sublimation (again in Figure 4, there are large counts for rain evaporation, snow sublimation, and the Bergeron process). These eigenvalues are again positive mainly when the affected particles are small and few in number, as such particles rapidly evaporate when out-of-cloud. In such cases, the particle mass rapidly drops to zero, and the temporal resolution should not be relevant to the final state of the system.
Turning back to the negative eigenvalues, we wanted to examine the processes associated with the smallest timescales. As seen in Figure 6, these processes are as follows: 1. Accretion of cloud water by rain. 2. Rain evaporation. 3. Rain self-collection. 4. Snow sublimation. 5. Snow self-collection. 6. Vapor/Ice transfer.
Other processes that were active with relatively short timescales include the Bergeron process and heterogeneous rain freezing, but we believe that the process rates are not being accurately calculated by our driver, and so we are less concerned about the apparent short timescales involved. The Bergeron process rate calculation is believed to be one of the least accurate parameterizations in MG2, primarily because it is poorly suited to the coarse spatial resolution used for GCMs, sometimes even resulting in a rate that is orders of magnitude too large (Tan & Storelvmo, 2016;Zhang et al., 2019). Reducing the time integration error is simply not a priority until this issue is addressed. As for the heterogeneous rain freezing, it is the process least Figure 5. Diagram of accretion regimes. In case (a), the cloud mass is still much larger than the rain mass. As the rain accretes more liquid mass, it becomes more effective at accreting additional mass, so accretion experiences positive feedback. In case (c), the rain has already accreted most of the cloud, and further depletion of cloud mass slows the accretion rate, producing negative feedback. In case (b), cloud and rain mass are comparable, and these effects are in balance, producing only weak positive or negative feedback. likely to behave the same way in our driver as it does in E3SM, because we have disabled instantaneous freezing of rain, and the heterogeneous process should be expected to be more active to compensate.
The particular processes that we examined in more detail were rain evaporation, self-collection, and accretion of cloud water. Partially this was because the purely liquid processes were easier to examine in isolation from the other physics, and partly this was because in practice the vapor/ice transfer was usually associated with small timescales only when the effect of this process was small anyway, so the error due to finite time resolution was small. In particular, sublimation of a small ice mass in relatively dry air can be quite rapid, but leads to the same behavior as given by the limited behavior, namely, that the ice sublimates completely and rapidly.
We can also note that almost no eigenvalues are associated with the external processes, which is what we expect since these process' rates do not directly depend on the MG2 inputs and only interact with MG2's physics due to being affected by the conservation limiters. In our data set, we found only three eigenvalues associated with contact freezing, a few thousand associated with immersion freezing, and none with nucleation deposition or droplet activation. Figure 7 shows the degree to which different processes were associated with the same eigenvalues (and hence timescales). Looking at timescales associated primarily with the processes listed above, we found the following: 1. Accretion of cloud water by rain is associated with autoconversion. 2. Rain self-collection is strongly associated with rain evaporation. 3. The degree of association between ice processes is quite complex. In particular, vapor/ice transfer is associated with several processes, possibly because it is mildly active in a very large number of grid cells.
As Figure 6 shows, the timescales associated with both liquid and ice autoconversion were relatively long. We therefore expected that the regime in which, say, liquid autoconversion and accretion interact most heavily would be different from the regime in which accretion by rain is associated with short timescales, and therefore we hypothesized that autoconversion should not require a short time step to adequately resolve, despite its association with accretion.
On the other hand, rain self-collection and evaporation were both primarily associated with the same short timescales, and so it appeared less likely that we could resolve rain-related processes without using a relatively short timescale for both. Figure 7. Association between pairs of processes in MG2. Each row shows the average association index (C) for eigenvalues associated with a given primary process. That is, each row of the table shows the average value of a row of C for eigenvalues primarily associated with the process shown on the left axis. By our definition of primary association, this means that all elements on the diagonal have values of at least 0.5, so each entry on the diagonal has had its value reduced by 0.5 to fit in the same color range.

Decomposition by Weather Regime
Microphysics operates differently in different meteorological conditions. In particular, only a few microphysical processes are typically active at any particular point in space and time. Thus, breaking cases down by weather regime can simplify the task of understanding model behavior. Such decomposition is also important because processes are likely to have different timescales depending on the weather regime they're in. For instance, the rate of accretion is fairly steady when rain and cloud mass are similar, but the accretion rate decreases rapidly when cloud becomes depleted, and the associated timescale is therefore much smaller in the latter case.

Methodology
We accomplished this decomposition by normalizing all process rates by typical values (so all processes were given approximately equal weight) and using a simple k-means algorithm to cluster grid cells based on process rates, so that we could treat the clusters as separate regimes.
We were interested in finding clusters containing qualitatively different, common types of behavior to examine in further detail, not in an objective categorization of all types of grid cells in the model. For this purpose we were content to use a degree of hand-tuning to determine both the total number of clusters and the scaling of process rates used in the clustering algorithm (for instance, reducing the effect of the Bergeron process). We iterated with changes to the cluster number and scaling until we found 10 clusters that had reasonably distinct process rates from one another.
Once this was accomplished, we could then look at the eigenvalues associated with the Jacobian for each cluster, in order to determine which regimes were likely to have short timescale behavior based on process rates alone. Table 2 gives a brief summary of the most active processes in each of these clusters, specifically those with mean rate above 10 −5 g/kg/s (or 0.298/kg/s for precipitation self-collection), as well as the number of grid points in each cluster. Note that these rates are 100 times greater than those used in section 3.2, since we are describing those processes that dominate the physics for each cluster, rather than those that simply have some appreciable activity. Note that Cluster 0, which consists of generally clear-sky grid cells where no processes are particularly active, comprises the majority of our data set.

Results
While the Bergeron process for cloud ice and heterogeneous rain freezing may seem quite large in many clusters, these should be disregarded for two reasons. First, as noted in section 4.2, these two processes are not likely to be well represented in our driver. For the Bergeron process, this is due to the low accuracy of Het. Rain Frz.: 3.63 × 10 −4 g/kg/s 4,579 Berg. (Cloud): 2.73 × 10 −5 g/kg/s Liq. Accr. Snow: 1.37 × 10 −5 g/kg/s Ice. Auto.: 1.00 × 10 −5 g/kg/s Rain Self-col.: 1.00/kg/s Snow Self-col.: 0.535/kg/s 9 Rain over a broad area. Rain Evap.: 3.48 × 10 −4 g/kg/s 29,040 Rain Self-col.: 44.5/kg/s MG2's representation of this process at the spatial resolutions used by GCMs. We believe that there will be little benefit to improving time integration of the Bergeron process until this more fundamental concern is addressed. As for the heterogeneous rain freezing, this process is likely affected by our stand-alone driver deviating from the behavior of MG2 within E3SM, in particular, by the fact that the instantaneous rain freezing has been disabled. Second, our scaling deemphasized these processes, so they had the less impact on the clustering algorithm itself than the other MG2 processes.
Note that this provides another, simpler way to look at associations between processes. Looking again at liquid autoconversion and accretion, we see that these two processes are active together in Clusters 3 and 6, but that the rate of accretion is much larger in Cluster 6. We can then look at the eigenvalues of the Jacobian based on cluster, which is shown in Figure 8. Unsurprisingly, the negative eigenvalues are much larger in Cluster 6, implying that the short timescales are associated with the heavy accretion present in cells with a relatively large in-cloud rain mass.
Similarly, we can use other clusters to examine other sets of short timescale processes. Cluster 1 would the best test case for short timescales associated with the Bergeron process and ice deposition (if the Bergeron parameterization was made more accurate). Rain evaporation and self-collection are both active and associated with short timescales in Clusters 7 and 9.

Impact of Shorter Time Steps
In the previous sections we identified processes and combinations of processes that evolve much more rapidly than the default model step. In this section we test whether accurately resolving those processes has a big impact on model behavior.

Methodology
In this section, we focus on two clusters, each associated with a set of processes that we believed would change considerably if substepped so as to resolve their behavior: 1. Grid cells with large rain self-collection and evaporation rates, corresponding to heavy out-of-cloud rain (Cluster 9). 2. Grid cells with large accretion rates and moderate autoconversion rates, corresponding to heavy in-cloud rain. (Cluster 6, filtered to remove grid cells with any ice process rate or rain evaporation rate above 10 −7 g/kg/s).
The subsampling of Cluster 6 was necessary to remove cases where accretion and autoconversion were not the only active processes. This was necessary for two reasons. First, the clustering algorithm alone did not perfectly separate those grid cells with only rain production from grid cells that combined rain production with evaporation/sublimation of precipitation falling from above. Second, in some grid cells, especially at large time step sizes, the cloud liquid is turned completely to rain. If all cloud mass is removed from a grid cell, MG2 no longer considers the fraction of the grid cell that was occupied by the cloud to be saturated, and rain evaporation turns on. Since we were interested in the effect of substepping on autoconversion and accretion specifically, we removed grid cells where this occurred from this portion of the study. It should be noted that rain self-collection is also active in Cluster 6. However, autoconversion and accretion within MG2 are not functions of the rain number, so it is not necessary to account for rain self-collection to accurately represent these processes. For each of these cases, we produced modified versions of the MG2 standalone driver from section 3.1.1, which allowed these processes to be run at a smaller time step using the forward Euler method. For instance, for heavy in-cloud rain, autoconversion and accretion were substepped inside a nested series of loops, so that it was possible to adjust the time step of each independently of MG2 as a whole. If both processes were run at a finer time step, it was also possible to independently adjust the coupling frequency. By adjusting these time steps, it was possible to determine which processes needed to be better resolved to improve MG2's accuracy, and which were less relevant.
In order to assess the accuracy of these substepped simulations, it was necessary to decide upon both a measure of error and a "converged" result for comparison. To measure the error in each grid cell, we used the total water mass difference defined above (D w ), which for Cluster 9 was effectively identical to the evaporation rate due to the absence of other influences (recall that D w only accounts for the error in water mass, not hydrometeor number). The converged result was chosen by using MG2 run at a very fine time step of 75/512 s, which is 1/2,048 of the normal MG2 time step of 300 s. Using a power of two in this way is convenient for convergence studies, where the time step can be progressively halved to provide results that are closer and closer to converged.
Note that for the forward Euler method, the time integration error should be proportional to the time step (first-order convergence). However, this is typically expected to hold only for short enough time steps, and in particular there is no guarantee that this will hold when the method is kept stable by limiters rather than because the method is absolutely stable at a given time step. Generally, there is little point in reducing the time step size if doing so does not much improve the accuracy, so when the time step for a process is reduced, we generally hope to reach the regime where first-order convergence holds, as a bare minimum.

Results
We hypothesized that for the grid cells with large rain evaporation, there would be a noticeable benefit from substepping rain self-collection as well, since our Jacobian-based analysis associates these two processes with the same short timescales.
Since these timescales are typically around 20 s, we also expect to see roughly first-order convergence for time steps shorter than this, but not for larger time step sizes, because for larger time steps this method (or rather, its linearization around a typical state) is not stable on this problem. The model therefore becomes increasingly dependent on limiters for longer time steps.
The actual evaporation/self-collection results are shown in Figure 9. For comparison, note that the mean evaporation rate in this cluster is about 1.82 × 10 −4 g/kg/s, so the mean error at a 300 s time step is an appreciable fraction of the evaporation rate itself. Substepping evaporation by itself is indeed much less effective than substepping both processes together. The self-collection of raindrops is a very fast process in MG2, and accounting for this considerably reduces the rate of evaporation. In particular, for grid points in Cluster 9, the self-collection causes rain drops to reach the maximum allowed size in less than 300 s, meaning that unless processes like rain evaporation are coupled at a smaller timescale, the particle size is effectively determined by the size limiter, not the details of the self-collection itself. This is one reason it is unsurprising that we see sublinear convergence at larger time step sizes.
Furthermore, the rate of convergence of the model becomes first order below a time step size of roughly 30-60 s, both when all of MG2 is substepped and when only these two processes are substepped together. In fact, it appears faster than first order when all of MG2 is substepped. This may be due to the fact that our "resolved" result is the 75/512 s result, and this causes some underestimation of the error in the leftmost few points on the graph. At very short time steps, the error levels off when substepping only evaporation and self-collection due to small contributions from other, nonsubstepped processes in the grid cell (e.g., small amounts of cloud are present, causing some accretion, though the average accretion rate in this cluster is orders of magnitude less than the evaporation rate).
For grid cells dominated by accretion and autoconversion, we expected the effect of substepping accretion to matter much more than that of autoconversion, since autoconversion is associated almost exclusively with timescales longer than the maximum MG2 step size of 300 s. Similar to the evaporation case, we expected that first-order convergence would occur only for timescales shorter than about 5 s. Figure 10 shows the results of this substepping of accretion and autoconversion in Cluster 6. Accretion rates are 1.13×10 −4 g/kg/s, so the error is again comparable to the accretion rate itself. Substepping autoconversion alone is ineffective for reducing the error. Substepping the accretion, however, can cut the errors by more than an order of magnitude, while substepping autoconversion as well improves the error by only another factor of 2.
However, we find, somewhat surprisingly, a near-perfect first-order convergence when substepping these processes, even out to 300 s. Note, again, the much smaller timescales in Cluster 6 shown in Figure 8, which would seem to suggest that this time step would produce instability, not linear convergence.
We came to suspect that by subsampling Cluster 6, we had removed the grid points with the shortest timescales. There are three possible reasons this could happen. First, the grid points with active ice processes could have been the grid points disproportionately associated with the shortest associated timescales. Second, the grid points with no ice processes but some evaporation at the beginning of the time step could be responsible, since some grid cells with out-of-cloud precipitation were included. Finally, some grid points Figure 11. Histogram of negative eigenvalues in Cluster 6. Blue bars correspond to eigenvalues at grid points excluded from the subsample due to having active ice processes (647 grid points). Green bars correspond to grid points excluded due to having active rain evaporation at the beginning of the time step (118 grid points). Yellow bars correspond to the grid points excluded due to developing rain evaporation over the course of a time step (3,796 grid points). Red bars correspond to the grid points included in the subsample (843 grid points). A black line is placed at (300 s) −1 for comparison with the MG2 time step.
could have no rain evaporation at the beginning of the time step, but develop rain evaporation over the course of a time step. In these grid points, at the beginning of the time step, all rain is considered to be inside the cloud, and no evaporation occurs. However, if the cloud is depleted to a concentration of less than 10 −3 g/kg, MG2 ignores the cloud area and assumes that all precipitation is outside the cloud, and therefore can evaporate. Figure 11 shows the results of our investigation. We find that our subsample did in fact include mostly points associated with long timescales, and excluded the majority of shorter timescales. Furthermore, the majority of points in Cluster 6, especially those associated with short timescales, were points that started with no active evaporation or ice processes, but which ended up depleting enough cloud to cause MG2's evaporation process to switch on.
Given this situation, we decided to use a second, broader subsample so as not to exclude points that develop rain evaporation. This corresponds to the yellow and red bars in Figure 11. For this broader subsample, substepping accretion and autoconversion alone would probably not be sufficient to adequately represent the rain mass overall, since rain evaporation also becomes an important process. However, we might hope that substepping accretion and autoconversion together would at least reduce error in the rain production rate, especially since most evaporation should occur after the cloud is depleted, when the accretion should be slowing down.
Unfortunately, this is not the case, as Figure 12 shows. There is significant error in the rain production rate unless all rain-related processes are substepped together. We note that first-order convergence happens below a time step of about 20 s. Figure 11 suggests that we might need an even smaller time step, but in practice, this does not seem to make much difference. This may be because these eigenvalues were calculated based on the conditions at the beginning of the time step. After some time has elapsed, the evaporation is turned on, while the self-collection is likely to cause the particles to reach maximum size, at which point the self-collection effectively turns off. This represents an important feature of MG2 not captured by our linear approach; if the time step is reduced, then different processes may activate at different times, and the timescales associated with active processes may also change over time.
It is worth noting that all results in this section are broadly consistent with the prior literature in that they suggest a time step size of 60 s or below is necessary to reduce the error to within a few percent of the total process rate. Notably, in the case of rain evaporation, the error hardly improves until the time step is on the order of tens of seconds or lower.

Discussion
Our analysis of MG2, based on examining the eigenvalues of a numerically derived Jacobian, shows that there are many situations where MG2 would appear to be unstable when using E3SM's forward Euler method at a 300 s time step. It is stable in practice due to nonlinearities in MG2, and especially due to the presence of limiters, but we would expect these limiters to reduce the accuracy with which E3SM solves the equations that MG2 is intended to implement.  This leads us to three concerns. First, is the time discretization error first-order, so that we can readily trade off increased computational cost for a proportional decrease in error? Or are there key unresolved timescales present in the system, so that the solutions we find at a 300 s time step are qualitatively different from the converged results? For instance, if certain process rates are typically constrained by limiters at long time steps, then modest improvements to the process rates will have little effect on the result, since both the original and "improved" process rates will be overwritten with the limited value. There will therefore be little benefit to reducing the time step until we reach the point where the model runs stably without these limiters.
Our results suggest that whether or not one considers MG2 to be first order depends significantly on the regime and the set of active processes involved. The timescales involved in the production and growth of snow, for instance, seem to be generally much longer than the 300 s time step MG2 typically runs at. For processes involving rain or cloud ice, however, MG2 relies on limiters for stability, and in cases where rain evaporation is an important process, we have shown that the 300 s time step is an order of magnitude too large to achieve first-order convergence.
Our second concern is whether we can use a process-based analysis of MG2 to improve accuracy without substepping the entire microphysics. Our results again suggest that it is possible to considerably improve accuracy in this way in certain regimes, but that it takes some care to do so effectively. For instance, substepping MG2's rain evaporation at a much shorter time step, say 10 s, produces a large increase in model cost Figure 15. Comparison of zonally averaged rain mass in E3SM runs using MG2 at a 300 s (left) and 1 s (right) time step over 3 years, with the difference also plotted (center).
(it requires multiple calculations involving computationally expensive exponential functions), yet this only moderately reduces the error. But by also substepping rain self-collection (a less computationally intensive process) in the same loop, the error can be reduced by an order of magnitude. In general, it seems necessary to substep rain-related processes together, but it may be possible to do so without substepping the entire microphysics. We believe that sedimentation may also be relevant to these rain processes, and hope that a future study will investigate including rain evaporation, self-collection, and accretion as part of the sedimentation solver itself, rather than sequentially splitting sedimentation from the microphysical processes.
Third, we can broaden our view and ask whether it is necessary at this point to accept the increased computational cost in order to improve MG2's accuracy at all. To answer this question, we ran a simulation using E3SM v1 with MG2 substepped at a 1 s time step (compset F1850C5AV1C-04P2, grid ne30_ne30) for four simulated years. Figure 13 shows a Taylor diagram (Taylor, 2001) that compares the spatial variability of several key variables in this run to a control run with the MG2 default time step of 300 s. The differences here appear to be quite minor. Figure 14 shows the spatial distribution of total precipitation, which likewise is quite similar and shows no obvious systematic differences. We do note that the global mean stratiform precipitation from MG2 increases from 1.23 mm to 1.31 mm per day. However, the convective precipitation is reduced to a degree that almost exactly cancels this, from 1.90 mm to 1.83 mm per day.
We can see much larger differences if we look specifically at the vertical distribution of rain mass from MG2. Figure 15 shows the difference between the zonally averaged in-area rain mass, which shows a reduction of rain in the upper cloud and near the surface, as well as an increase in the lower cloud. This is due to a combination of reduced production of rain, reduced evaporation, and the effects of these differences on the rain fall speed. In particular, if a reduction in rain evaporation in the cloud base is matched by an equal reduction in rain production within the cloud (both from the stratiform and convective schemes), that suggests that the model is not particularly sensitive to the rain evaporation rate, since the precipitation is governed largely by top-down constraints rather than such specific details of the microphysics. However, the exact mechanism for this feedback is not immediately clear; one could also easily imagine an alternative outcome where a change in rain process rates affects cloud properties enough to cause a significant change in the shortwave cloud forcing. We also believe that the coupling of sedimentation to these processes should be further examined.

Conclusions
By numerically calculating the eigenspectrum of the Jacobian of the MG2 microphysics, we were able to associate a set of timescales to these microphysical processes, and found that MG2's equations are not always accurately modeled at the default E3SM time step. We were able to associate the short timescales to the coupled behavior of rain evaporation and self-collection for some grid cells, and demonstrated that in some regimes, decreasing the MG2 time step, particularly below 60 s, can lead to notable differences in rain-related processes. We have also demonstrated that rain evaporation and self-collection should not be treated independently in MG2's formulation, though accretion appears to be more independent.