Improving Time Step Convergence in an Atmosphere Model With Simplified Physics: Using Mathematical Rigor to Avoid Nonphysical Behavior in a Parameterization

Global atmospheric models seek to capture physical phenomena across a wide range of time and length scales. For this to be a feasible task, the physical processes with time or length scales below that of a computational time step or grid cell size are simplified as one or more parameterizations. Inadvertent oversimplification can violate constraints or destroy relationships in the original physical system and consequently lead to unexpected and physically invalid behavior. An example of such a problem has been investigated in the work of Wan et al. (2020, https://doi.org/10.1029/2019MS001982). This work addresses the issues at a more fundamental level by revisiting the parameterization derivation. A derivation of an unaveraged condensation rate in the unaveraged equations, sometimes referred to as subgrid equations, provides a clear description and more accurate quantification of the condensation/evaporation processes associated with cloud growth/decay, while avoiding simplifications used in earlier studies. A subgrid reconstruction (SGR) methodology is used to connect the unaveraged condensation rate with the grid cell averaged equations solved by the global model. Analyses of the SGR method and the numerical results provide insights into root causes of inconsistent discrete formulations and nonphysical behavior. It is also shown that the SGR methodology provides a flexible framework for addressing such inconsistencies. This work serves as a demonstration that when nonphysical behavior in a parameterization of subgrid variability is avoided through rigorous mathematical derivation, the resulting formulation can exhibit both better numerical convergence properties and significant impact on long‐term climate.


Introduction
The modeling of global atmospheric circulation is a difficult multiscale endeavor. For climate science, an effective atmospheric general circulation model (AGCM) must consider both small and large spatial and temporal scales: from the behavior of water droplets (millimeters and seconds) to behavior of clouds (meters and hours) and on to rainfall patterns that extend thousands of kilometers and progress on time scales of days to months. With current and upcoming computational systems, the computational grids in today's flagship AGCMs are on the resolution of a few tens to a hundred kilometers. For example, the current high-and low-resolution configurations of the AGCM in Energy Exascale Earth System Model (EAMv1, Rasch et al., 2019) have 1°and 1/4°horizontal resolution, respectively. Even when the model reaches its next target of 3 km resolution in 5 to 10 years (E3SM Project, 2020), the grid spacing will still be too coarse to directly resolve certain physical processes, including cloud microphysics (Pruppacher & Klett, 2000, Figures 2-4 and 2-5) and turbulence ( Stull, 1988, Figure 1.15). For such physical processes, parameterizations must be used to describe the impact on larger scales.
Cloud processes are one class of atmospheric phenomena for which parameterizations are necessary in current global models. Clouds have essential impacts on both the energy and hydrological cycles in the atmosphere. Warm (liquid) clouds are formed when cloud droplets are created from the condensation of water vapor. Water vapor condensation occurs due to the cooling and/or moistening of the atmosphere from upward motion, radiation, and advection. The condensation associated with nonconvective clouds is commonly referred to as large-scale condensation in global models (Park et al., 2014;Rasch & Kristjánsson, 1998;Sundqvist, 1978;Tiedtke, 1993;Wilson et al., 2008).
Deriving a parameterization for a given physical process is not a trivial task (McFarlane, 2011). As described by Simon and Behrens (2019), parameterizations derived in a heuristic or intuitive manner can lead to incorrect and unintended physical behaviors. Park et al. (2014) provide an example for atmospheric parameterizations, discussing the inconsistency between diagnosed nonconvective cloud fraction and prognosed cloud condensate amount. Numerical remedies such as clipping limiters, which replace values that exceed predetermined bounds by a predetermined acceptable value (e.g., negative values replaced by 0), can mask the existence of problematic behaviors and leave the causes of those behaviors unaddressed. Zhang et al. (2018) showed that clipping limiters can result in spurious sources and sinks of atmospheric water with nonnegligible magnitudes. Wan et al. (2013) showed that the underlying numerical issues that triggered clipping limiters for the sulfuric acid gas concentration in a climate model led to significant impacts on the simulated aerosol life cycle. In the companion paper to this work, Wan et al. (2020) showed that a large-scale condensation parameterization derived from assumptions on grid cell average quantities following Rasch and Kristjánsson (1998) and M. Zhang et al. (2003) leaves room for inconsistency between the simulated cloud fraction and liquid condensate amount. That inconsistency can result in large, nonphysical values of in-cloud liquid concentration. Wan et al. (2020) also showed that although the most egregious cases occur only in a small fraction of the grid cells, they can cause the measured time step convergence of the numerical solutions to degrade substantially from the expected value. Changes in the model formulation and time stepping method that address such nonphysical behaviors can substantially improve solution convergence in short (less than a day) simulations and also significantly change the statistics of the simulated long-term climate. Compared with Wan et al. (2020), who addressed the nonphysical behaviors by intuitively revising a closure assumption and improving the splitting of the time integration, this work addresses the issues at a more fundamental level by revisiting the parameterization derivation.
While it is common to construct a parameterization using assumptions on grid cell average quantities (Rasch & Kristjánsson, 1998;Sundqvist, 1978;Tiedtke, 1993;M. Zhang et al., 2003), an alternative approach is to first derive a formulation directly from a description of the physical phenomena that enforces the desired physical relationships and constraints, and then choose a numerical method to establish connections between this formulation and the grid cell average quantities the AGCM solves for. For example, both the reconstruct-evolve-average (LeVeque, 2002) and semi-Lagrangian multiscale finite element method (Simon & Behrens, 2019) approaches reconstruct spatial profiles within a grid cell from discrete grid cell average values. Such profiles, referred to as subgrid profiles, can be evolved using equations independent of grid cell size and then be mapped back to discrete grid cell average values. This type of approach, herein referred to as the subgrid reconstruction (SGR) methodology, can locally enforce conservation and constraints and also preserve asymptotic behavior.
This work applies the SGR methodology to the the large-scale condensation/evaporation processes discussed in Wan et al. (2020). Rigorous mathematical derivation provides an unaveraged version of the condensation rate that clearly describes the physical processes associated with condensation/evaporation in different scenarios. Combined with reconstruction and averaging to provide the grid cell averages needed by the AGCM, this derivation both avoids a simplification introduced by Rasch and Kristjánsson (1998) and later used by M. Zhang et al. (2003) and Wan et al. (2020) and inherently avoids the singularity problem that caused nonphysical behavior and a convergence problem in Wan et al. (2020). The flexibility of the SGR approach, specifically in the choice of subgrid profiles, can be used to address inconsistencies.
The SGR methodology has been used in the atmospheric modeling community (Carpenter et al., 1990;Kaas, 2008;Lauritzen et al., 2010;Malardel & Ricard, 2015;Zerroukat et al., 2002). Unlike previous work in the literature that often uses the SGR methodology for accurate simulation of fluid dynamics or tracer transport, the emphasis here is on avoiding pathology and ensuring proper convergence in a parameterization. Despite the fact that the rate of time step convergence has not been used as a standard metric in the evaluation of physics parameterizations, especially in the context of global simulations, achieving the expected convergence can help to build confidence that the simulations correspond to a correct implementation of the numerical methods. In addition, as is shown here and in Wan et al. (2020), poor convergence can be an indicator of nonphysical behavior, and addressing convergence issues can lead to changes in the long-term climate of an AGCM that are sizable and nonnegligible to atmospheric model developers. Finally, proper convergence rates allow global atmospheric models to achieve the expected accuracy gain from time step and grid refinement opportunities provided by upcoming increases in computational power. While an SGR method that requires values from many neighboring grid cells may be more computationally expensive than a parameterization that only needs the local grid cell average, the compact SGR method used in this study does not significantly contribute to the computational cost, verified by timing measurements.
The remainder of the paper is organized as follows. Section 2 provides an overview of the SGR methodology as used in this work. Section 3 derives the unaveraged condensation rate. Section 4 describes the reconstruction of subgrid profiles and the averaging that provides the grid cell average condensation rates. Several potential subgrid reconstruction profiles are discussed, each leading to its own variant of the grid cell average condensation rate. The temporal discretization of the averaged equations are presented in section 5. In section 6, the numerical results are presented and two evaluation criterion are discussed: solution self-convergence with respect to time step size and climatological effect. Guided by those criteria, one particular variant of the grid cell average condensation rate is found that more properly describes the underlying physics while having significant climatological impact. Section 7 summarizes the work and draws the conclusions.

Overview of Methodology
As in the companion paper by Wan et al. (2020), this work demonstrates the key methodology using a simplified global AGCM that considers only the fluid dynamics (dynamical core) and its interaction with the large-scale condensation of water vapor (or reversely, the evaporation of cloud liquid). Other physical processes typically included in an AGCM, such as radiation and cloud microphysics, are not present. The key effects of water vapor condensation include the conversion of water mass from gas phase to liquid phase and the release of latent heat. Evaporation of liquid water has the opposite effects. Following Rasch and Kristjánsson (1998) and M. Zhang et al. (2003), the time evolution equations for air temperature T(x, t) (unit: K), water vapor concentration q v (x, t) (unit: kg/kg), and cloud liquid concentration q l (x, t) (unit: kg/kg) at any spatial location x can be written as (1) where Q(x, t) (unit: kg/kg/s) is the condensation rate. The quantities A T (x, t) (unit: K/s), A v (x, t) (unit: kg/kg/s), and A l (x, t) (unit: kg/kg/s) denote the rate of changes in the atmospheric state that are caused by processes external to the phase change processes. Following the convention in EAM and its predecessors, the tracer concentrations are defined as the mass of tracer per unit mass of moist air (dry air plus water vapor), and the tendencies are defined accordingly.
Because of the practical need to solve equation set (1)  (2) where the overline notation denotes a grid cell average quantity. Wan et al. (2020) showed that if an expression is derived for QðtÞ directly from other grid cell-averaged quantities with overly simplified assumptions, unexpected or even physically invalid behavior can occur. This concern motivates the use of the unaveraged Equation 1 to first derive an expression for Q(x, t) that can then be averaged over a grid cell to obtain QðtÞ in (2). Furthermore, using (1) allows for the application of traditional numerical analyses. This paper presents such a derivation via the following steps:

Model Equations for (Unaveraged) Condensation
To derive an expression for Q(x, t), the equations in (1) are integrated from t to t + τ, where τ is an arbitrary positive number (a time interval), to give Note that the A * ½ · ðx; t; τÞ terms represent a time-averaged version of the external effects A [·] (x, s). In the spirit of the work by Sundqvist (1978), such external effects are cooling/warming and moistening/drying mechanisms in the atmosphere that trigger condensation/evaporation. As such, the A * ½ · ðx; t; τÞ terms can be viewed as forcing terms for a condensation parameterization.
While a mathematical description of water vapor condensation can be obtained by explicitly describing the mass transfer at the air-liquid interface of each water droplet, the spatial and temporal scales of such a model would not be suitable for a global atmosphere simulation conducted using typical spatial and temporal resolutions. The alternate that is typically utilized by global AGCMs is to assume that condensation happens immediately after the air is saturated and a new equilibrium is reached instantaneously (e.g., Golaz et al., 2002;Park et al., 2014;Sundqvist, 1978;Tompkins, 2002;Wilson et al., 2008). Here, the same assumption is made. Denote q sat (x, t) as the saturation vapor concentration. Under the instantaneous equilibrium assumption, q v (x, t) ≤ q sat (x, t). If q v (x, t) < q sat (x, t), it is assumed that x is cloud-free at time t in the sense that q l ðx; tÞ ¼ 0. If x is cloudy at time t (i.e., q l (x, t) > 0), it is assumed that q v ðx; tÞ ¼ q sat ðx; tÞ. These constraints on q v and q l , namely, q v ðx; tÞ ≤ q sat Tðx; tÞ; pðx; tÞ ð Þand q v ðx; tÞ < q sat Tðx; tÞ; pðx; tÞ ð Þ ⇒ q l ðx; tÞ ¼ 0; are used to determine the condensation rate Q(x, t) in (3)-(5).
Next, the atmospheric state T(x, t), q v (x, t), and q l (x, t)), as well as the forcing terms A [·] (x, t), are used to determine which of the four scenarios is happening at location x at time t (section 3.2) and what the condensation rate is for each scenario (section 3.1).

Scenario Condensation Rates 3.1.1. Scenarios PF and CA (Cloud-Free at t + τ)
Both the persistent cloud-free and cloud annihilation scenarios require that q v (x, t + τ) < q sat (x, t + τ) by definition. The liquid concentration constraint (9) gives that q l ðx; t þ τÞ ¼ 0 for these scenarios. Using this constraint in (5) gives The physical interpretation of Equation 10 is that, in order to create a cloud-free scenario at t + τ, evaporation is needed to remove all the cloud liquid present at t as well as any "new" liquid that will be brought to the location by other physical processes.

Scenarios PC and CF (Cloudy at t + τ)
Both the persistent cloudy and cloud formation scenarios require that q v ðx; t þ τÞ ¼ q sat ðx; t þ τÞ. For these two scenarios, an expression for q sat (x, t + τ) is needed in terms of quantities at time t.
Note that this Taylor series approximation requires that q sat is twice continuously differentiable in T and p and that pðx; t þ τÞ − pðx; tÞ ð Þ ð ∂q sat =∂pÞðx; tÞ is negligibly small (i.e., Oðτ 2 Þ). The same assumptions have been used in the simple-physics parameterization suite of Reed and Jablonowski (2012) and in the parameterization by M. Zhang et al. (2003) that was used in CAM versions 2-4. Assuming these conditions are satisfied, the condensation rate for the persistent cloudy and cloud formation scenarios is now derived from q v ðx; t þ τÞ ¼ q sat ðx; t þ τÞ, the Taylor series approximation, and the vapor concentration Equation 4: where Referring to the expression q sat (x, t) − q v (x, t) as the saturation deficit at time t and recalling that the scenarios under consideration are cloudy at t + τ, the physical interpretation of Equation 14 is that, if the air is already saturated at time t (scenario PC), then condensation will occur and compensate the external moistening and/or cooling effects in order to keep q v at saturation. If the air has a saturation deficit at time t (scenario CN), then external moistening and cooling that can overcome the saturation deficit will lead to condensation.

Scenario Determination 3.2.1. Scenarios PF and CF (Cloud-Free at t)
Consider now when q v (x, t) < q sat (x, t). Either the persistent cloud-free or cloud formation scenario could be occurring, so it is not immediately clear which of (10) or (13) to use for Q * (x, t, τ). If the persistent cloud-free scenario is considered, the condensation rate is (10) and it is expected that q v (x, t + τ) < q sat (x, t + τ): or, in other words, , the persistent cloud-free scenario is not possible, and the cloud formation scenario must be occurring. Additionally, note that if the cloud formation scenario is considered when , the cloud formation scenario is not possible, and the persistent cloud-free scenario must be occurring. Thus, whether or not D( The physical interpretation of the above derived scenario determination criterion is the following: because the air is cloud-free at time t (meaning q l (x, t) = 0), the quantity D * (x, t) reduces to the saturation deficit q sat (x, t) − q v (x, t). The quantity D(x, t, τ), on the other hand, is the cooling and moistening effects occurring during the time interval τ (including the import of cloud liquid that could potentially evaporate and also lead to moistening and cooling). The air becomes saturated and cloudy if the moistening and cooling effects can overcome the saturation deficit (D(x, t, τ) > D * (x, t)) or stays unsaturated (and cloud-free) otherwise (D(x, t, τ) < D * (x, t)).

Scenarios PC and CA (Cloudy at t)
Consider now when q v ðx; tÞ ¼ q sat ðx; tÞ. Either the persistent cloudy or cloud annihilation scenario could be occurring, so it is again not immediately clear which of (10) or (13) to use. If the cloud annihilation scenario is considered, the condensation rate is (10), and it is expected that q v (x, t + τ) < q sat (x, t + τ). It has already been shown that this condition requires D(x, t, τ) < D * (x, t), so that q v ðx; tÞ ¼ q sat ðx; tÞ and D(x, t, τ) ≥ D * (x, t) means the persistent cloudy scenario is occurring. If the persistent cloudy scenario is considered 10.1029/2019MS001974

Journal of Advances in Modeling Earth Systems
, then q l (x, t + τ) will again be less than 0. Thus, whether or not D(x, t, τ) < D * (x, t) also uniquely determines which scenario is occurring when q v ðx; tÞ ¼ q sat ðx; tÞ.
From a physical perspective, the PC and CA scenarios are cloudy at time t and hence saturated, giving D * ðx; tÞ ¼ −γðx; tÞq l ðx; tÞ. Whether the air can become cloud-free at t + τ depends on whether the external drying and/or warming effects can consume all the existing liquid (D(x, t, τ) < D * (x, t)).

Summary
The description of the condensation process happening at location x between t and t + τ is summarized in Table 1. Note that if the differential Equation 1 are preferred over the integral Equations 3-5, then Q(x, t) can be obtained by taking lim τ→0 Q * ðx; t; τÞ (see Appendix A).

Reconstruct and Average
As mentioned in section 2, the model Equations 3-5 can be averaged over a computational grid cell Ω to obtain grid cell average model equations (recall the overline notation indicates a grid cell average): One of the key features that distinguish this work from Wan et al. (2020), M. Zhang et al. (2003) and Rasch and Kristjánsson (1998) is that the evaluation of Q * ðt; τÞ is done by direct averaging of Q * (x, t, τ), defined in Table 1, across all the subgrid locations in the grid cell (i.e., x ∈ Ω). The various expressions for Q * (x, t, τ) require other subgrid values (e.g., q v (x, t), q l (x, t)). As mentioned in section 2, the AGCMs expect to evolve grid cell average quantities and, therefore, do not provide the required subgrid quantities. As such, spatial profiles for the subgrid quantities must be reconstructed from their corresponding grid cell average counterparts: pðtÞ; T ðtÞ; q v ðtÞ; q l ðtÞ; pðx; tÞ; Tðx; tÞ; q v ðx; tÞ; q l ðx; tÞ Once the subgrid quantities are available to determine Q * (x, t, τ), the grid cell average condensation rate Q * ðt; τÞ is obtained by averaging Q * (x, t, τ). In this paper, this reconstruct and average process is referred to as subgrid reconstruction.

Reconstruct
The simplest approach to subgrid reconstruction is to assume quantities are uniformly distributed throughout the grid cell; however, such a simplification to certain quantities can lead to unacceptably large errors at current grid cell resolutions. An example familiar to atmospheric scientists would be the sudden formation or annihilation of a 100 km by 100 km cloud. How to represent subgrid variations is a research topic by itself  (17))

Journal of Advances in Modeling Earth Systems
(see, e.g., Griffin & Larson, 2016;Larson et al., 2001Larson et al., , 2002K. Zhang et al., 2016). It is even less obvious how to represent quantities that do not have ample real-world measurements, for example the tendency termsA * ½ · in our study. In the following, some approximations from the earlier work of Rasch and Kristjánsson (1998) and M. Zhang et al. (2003) are adopted, but we use a more general formulation for q v , q l , A * l to allow for the exploration of more options.
The assumption of uniform subgrid distribution is applied to temperature, the temperature tendency, and pressure: Tðx; tÞ ¼ T ðtÞ, A * T ðx; t; τÞ ¼ A * T ðt; τÞ, and pðx; tÞ ¼ pðtÞ, following M. Zhang et al. (2003). It follows that the saturation water vapor concentration is also uniform across the grid cell: q sat ðx; tÞ ¼ q sat T ðtÞ; pðtÞ À Á , denoted simply as q sat ðtÞ for notation brevity. Similarly, The concept of cloud fraction is used to take into account the subgrid variation in cloud cover, as is commonly done in global AGCMs. Let f(t) denote the fraction of a grid cell that contains clouds at time t. A simplified version of the cloud fraction scheme used in CAM4 (Neale et al., 2010(Neale et al., , 2013 is used here for f(t): where RH ðtÞ: ¼ q v ðtÞ=q sat ðtÞ is the grid cell average relative humidity with RH * ¼ 0:8.
While x represents an arbitrary location in a three-dimensional geometry, one-dimensional cells are used for reconstruction. This is justified similar to how cloud fraction is used despite only describing the abundance of cloudy locations within the area of interest. Cloud fraction provides no information on the geographical coordinates of the subgrid locations within each two-dimensional horizontal grid cell nor on the variation of cloud amount in the vertical direction within each layer of the computational mesh.
Consider Ω ¼ ½0; L that is cloudy for x ∈ [0, c(t)) and cloud-free for x ∈ (c(t), L], where cðtÞ ¼ f ðtÞL. The cloudy region of the grid cell must be saturated: q v ðx; tÞ ¼ q sat ðtÞ for x ∈ [0, c(t)). Consider the vapor concentration to be a polynomial of degree p q v in the cloud-free region: Note that this is the unique piecewise-defined polynomial q v ðx; tÞ ¼ q sat ðtÞ for x < c(t) and q v ðx; tÞ ¼ aðtÞþ The cloud-free region of the grid cell must be free of cloud liquid: q l ðx; Consider the liquid concentration to be a polynomial of degree p q l in the cloudy region: Similarly to the vapor concentration polynomial, this is the unique piecewise-defined polynomial that satisfies q l ðtÞ ¼ ð1=LÞ The tendencies of vapor and liquid concentration are not governed by the constraints of (8) and (9). As such, the vapor concentration tendency is chosen to be spatially constant across the grid cell. The liquid concentration tendency is treated differently: while a uniform subgrid tendency seems harmless when the liquid amount increases in a grid cell (the average tendency is positive), a negative tendency needs to be 10.1029/2019MS001974

Journal of Advances in Modeling Earth Systems
attributed to the cloudy region in order to remain physical. As such, the cloud liquid tendency is considered a polynomial that is 0 at x ¼ L to focus the effect of the tendency on the cloudy region. Thus, the vapor and liquid concentration tendency polynomials are All the polynomial profiles are shown in Figure 1 for and quartic (p ½ · ¼ 4) polynomials.

Average
The grid cell average condensation rate is now obtained by averaging quantities from the one-dimensional reconstruction. Ideally, the determination criterion specified in Table 1 would be used to identify which of the persistent cloudy, persistent cloud-free, cloud formation, or cloud annihilation scenarios is occurring at each location x. The corresponding expressions of Q * (x, t, τ) would then be averaged via integrating from The challenge of this approach is that it would require identifying regions where q v ðx; tÞ < q sat ðtÞ and where D( Assuming all the roots of D(x, t, τ) − D * (x, t) can be found, the complexity of then considering all possible combinations of sign changes of D(x, t, τ) − D * (x, t) across roots and across x = c(t) makes it very unlikely that a closed-form expression for Q * ðt; τÞ can be obtained for use in an AGCM. For these reasons, the one-dimensional grid cell is instead partitioned into regions using c(t) = f(t)L and c(t + τ) = f(t + τ)L with the expression of f empirically chosen based on physical arguments (21). Note that the f(t + τ) given by this empirical expression might contradict the cloud fraction corresponding to D(x, t, τ) − D * (x, t) in the subgrid spatial profiles. Further discussion on this can be found in sections 6.2 and 6.3.
The cases where the cloud fraction of a grid cell decreases or increases from t to t + τ are considered separately, as these two cases can be further decomposed into different collections of scenarios defined in section 3. First, consider cases of cloud fraction decrease, when f(t + τ) < f(t), so that the grid cell has the persistent cloudy scenario for x ∈ [0, c(t + τ)), the cloud annihilation scenario for x ∈ (c(t + τ), c(t)), and the persistent cloud-free scenario for x ∈ (c(t), L]. The grid cell average condensation rate is now obtained by combining spatial integrals under the three scenarios: which gives where c A * l t; τ; p Al Term A represents the condensation/evaporation needed to keep the specific humidity in the cloudy region at saturation. Term B represents the condensation/evaporation needed to keep the cloud-free region free of cloud liquid, with c A * l t; τ; p Al À Á representing the average of the liquid concentration tendency profile in the cloud-free region. Note that this average decreases as p Al increases, which is consistent with higher-order polynomial profiles concentrating more cloud liquid tendency outside of the cloud-free region (into the cloudy region). Term C represents the evaporation needed to remove cloud liquid from the cloud annihilation region, with b q l t; p q l representing the average of the liquid concentration profile in the annihilation region. Note also that because f (t + τ) < f (t) implies that 1 − f (t + τ)/f (t) < 1 in (27) , this average b q l t; p q l decreases as p q l increases, which is consistent with higher-order polynomial profiles concentrating more liquid concentration outside of the annihilation region (into the persistent cloudy region). To compare (26) to the grid cell average condensation rate studied in Wan et al. (2020), consider using piecewise-defined constant polynomials for liquid concentration and concentration tendency p Al ¼ p q l ¼ 0 . In the limit as τ → 0, both grid cell average condensation rates are identical. It is worth noting that this expression of grid cell average condensation rate inherently avoids the singularity problem associated with the C term in the baseline model in Wan et al. (2020): the expression (26) applies only to cloud fraction decrease situations, in which case f (t) in the denominators in (27) will be nonzero. Now consider cases of cloud fraction increase, when f(t + τ) > f(t). The grid cell average condensation rate in the one-dimensional cell is which gives 10.1029/2019MS001974

Journal of Advances in Modeling Earth Systems
Terms A and B retain the same interpretation as in the cloud fraction decrease cases discussed earlier. Term C now represents the amount of condensation compensating for the portion of external moistening and/or cooling effect that would otherwise lead to supersaturation in the formation region. The value c sat t; τ; p q v is the average saturation deficit over the cloud formation region. This average decreases as p q v increases, which is consistent with higher-order polynomial profiles concentrating unsaturated water vapor in the formation region. To compare with Wan et al. (2020), consider piecewise-defined constant polynomials for vapor concentration and liquid concentration tendency p Al ¼ p q v ¼ 0 . In the limit as τ → 0, the two grid cell average condensation rates have identical A and B terms. The C term here has a clear physical meaning, namely the moistening and cooling that overcome saturation deficit leads to condensation. In contrast, the C terms in Rasch and Kristjánsson (1998), M. Zhang et al. (2003), and Wan et al. (2020), assume newly formed clouds have the same in-cloud liquid concentration as the existing clouds, which is clearly a simplification with room for improvement.

Temporal Discretization
Denote the computational time step as Δt. For a given discrete time t n , denote the next discrete time as t n þ 1 ¼ t n þ Δt. The A * ½ · terms in our simplified model (1) contain only the fluid dynamics and resolved-scale transport, so the approximations are provided by the dynamical core (in a sequentially split fashion) and denoted A * T n , A * v n , and A * l n . Denote Q * n as the first-order approximation to Q * ðt n ; ΔtÞ where all OðΔtÞ terms are neglected:

Journal of Advances in Modeling Earth Systems
The quantities T n þ 1 , q v n þ 1 , and q l n þ 1 are defined as The superscript terms T m , q v m , and q l m , with m being either n or n + 1, are approximations to T ðt m Þ, q v ðt m Þ, and q l ðt m Þ , respectively. The approximations are first order due both to the sequential splitting and to the dropping of OðΔtÞ terms. Note that the grid cell average vapor concentration q sat ðt n Þ and its derivative ∂q sat ∂T ðt n Þ are approximated as q sat n : ¼ q sat T n ; p n À Á and ∂q sat ∂T n : ¼ ∂q sat ∂T T n ; p n À Á , respectively, with the specific expression for q sat being the same as in Wan et al. (2020).

Approximating the Cloud Fraction at t n and t n + 1
Recall the cloud fraction scheme (21) is a function of the grid cell average relative humidity RH ðtÞ: ¼ q v ðtÞ =q sat ðtÞ. As such, the quantity f(t n ) is easily approximated as f n : ¼ f q v n =q sat n ð Þ . The cloud fraction at t n + 1 is more difficult to approximate. It is worth noting that df/dt is a function of dRH =dt, which depends on the grid cell average condensation rate. Thus, the value for f(t n + 1 ) cannot be known a priori when computing Q * n in (31). As such, a second-order extrapolation scheme is used to predict f(t n + 1 ) from quantities known at t n : 2. Compute an initial Q * n , denoted f Q * n , using e f n þ 1 . 3. Compute temporary updates for temperature and vapor concentration.
4. Compute a corrected e f n þ 1 using f e q v =q sat ð e T ; p n Þ in (21).

5.
If corrected e f n þ 1 equals 0 or 1, then compute corrected f Q * n using corrected e f n þ 1 .

Set
Note that the additional error from this correction step can be shown to be OðΔtÞ and, therefore, is better than ad hoc fixes that might cause large Δt-independent errors.

Implementation
The condensation parameterization is implemented in a revised version of the Energy Exascale Earth System Model (E3SM) code (E3SM Project, E3SM Project), using the hydrostatic spectral-element dynamical core for the resolved-scale fluid dynamics and tracer transport (Dennis et al., 2012;Taylor & Fournier, 2010). The initial state is defined as follows: 10.1029/2019MS001974

Journal of Advances in Modeling Earth Systems
• Set T 0 and q v 0 to values taken from a multiyear climate simulation conducted with the E3SM atmosphere model (EAM, Rasch et al., 2019), so that the initial temperature and vapor concentration values are a reasonable representation of realistic meteorological conditions. • Set q l 0 ¼ 0 to account for a different cloud fraction scheme used in EAM than here. • The remaining prognostic variables needed by the spectral-element dynamical core (wind velocity, surface pressure, etc.) are set in a similar manner as T 0 and q v 0 .

Numerical Results
The remainder of this paper focuses on evaluating the convergence and climatological effect of the subgrid reconstruction (SGR) grid cell average condensation rate derived in sections 3-5. Recall that the expression for the grid cell average condensation rate Q * in (31) has three parameters: p q v , p q l , and p Al . These parameters respectively determine the polynomials used to reconstruct subgrid profiles for q v , q l , and A l from the corresponding grid cell averages (q v , q l , A l ). For a given combination of p q v , p q l , and p Al , the corresponding grid cell average condensation rate is denoted SGR-P p q v ; p q l ; p Al h i . Wan et al. (2020) found that abnormal self-convergence (Knupp & Salari, 2002) behavior (degraded convergence rate) was associated with physically invalid behavior, specifically cloud-free grid cells with infinitely large in-cloud liquid concentration. Furthermore, Wan et al. (2020) found that changes in the parameterization and time discretization that helped improve convergence in short simulations also had climatological effects: affecting total cloud cover, longwave cloud radiative effect, and shortwave radiative effect. Both of these findings motivate the evaluation of solution features both in short diagnostic runs and in multiyear climate simulations with various condensation rates SGR-P[p q v , p q l , p Al ]. The results from the piecewiseconstant reconstruction (SGR-P[0, 0, 0]), which uses piecewise-constant quantities following M. Zhang et al. (2003) and Wan et al. (2020) motivates further investigation in the piecewise-linear reconstruction (SGR-P[1, 1, 1]) and piecewise-linear-quartic reconstruction (SGR-P[1, 1, 4]).

Experimental Design
For the self-convergence testing, the result from a simulation using time step size Δt ¼ 1 s is used as a reference solution, denoted with a superscript asterisk (e.g., T * ðtÞ is the reference temperature solution). With a subscript denoting the value in a certain grid cell (e.g., T i ðtÞ denoting the temperature value in the ith grid cell), the error is then computed with the same weighted root-mean squared (WRMS) error norm used in Wan et al. (2020): Δt is the time step used in obtaining T ðtÞ, dA is the horizontal area of the grid cell, and dp is the grid cell average pressure change across the grid cell. The values of Δt used are 4,8,15,30,75,120,450, and 1,800 s. The convergence rate itself is computed using a method similar to Wan et al. (2020). A linear, least squares fit is performed on fðlogðΔt k Þ; logðe k ðt; ΔtÞÞg for Δt k ∈ [8 s, 120 s]. This range is chosen a posteriori so that the time steps {Δt k } are in the intermediate asymptotic regime where the error exhibits a power law dependence on Δt. Note that if an analytic solution was available instead of the Δt ¼ 1 s reference, the lower bound on this range would be lower. The slope of the resulting linear, least squares fit is taken to be the convergence rate.
To study the climatological effects of using various SGR condensation rates, the condensation rates were implemented in the CAM4 physics parameterization suite which has the scheme of M. Zhang et al. (2003) as the default. 10-year climate simulations were conduced at approximately 1°horizontal resolution forced by climatological sea surface temperature. The zonal and annual mean total cloud cover and cloud radiative effects are shown below. A detailed description of the CAM4 physics suite can be found in Neale et al. (2010).

10.1029/2019MS001974
Journal of Advances in Modeling Earth Systems

Piecewise-Constant Reconstruction Results (SGR-P[0, 0, 0])
A self-convergence test is conducted for the piecewise-constant subgrid reconstruction for 1, 2, and 3 hr simulations. The resulting WRMS error (38) in temperature and the self-convergence rates are shown in the right panel of Figure 2. The dashed black line is a first-order reference line. The remaining lines show the linear least squares fit described earlier in this section, with the solid portions indicating time steps used for the fit. For comparison, the left panel of Figure 2 shows the results from the parameterization described in section 3.3 of Wan et al. (2020), herein referred to as the "baseline" parameterization. The measured convergence rate from SGR-P[0, 0, 0] is 0.9 after 1 hr and 0.8 after 2 or 3 hr, showing significant improvements over the convergence rates of 0.4 from the baseline simulations. Recall that a substantial difference between the baseline parameterization and the SGR approach is that the baseline parameterization uses the same grid cell average condensation rate expression regardless of whether the cloud fraction is increasing or decreasing. The SGR grid cell average condensation rate (31) has a term whose definition depends on whether f(t n + 1 ) ≥ f(t n ). Thus, discerning between cloud formation (f n + 1 > f n ) and cloud annihilation (f n + 1 < f n ) in (31) fundamentally avoids the problem of having zero cloud fraction in the denominator and hence shows better convergence over the baseline parameterization (see Wan et al. (2020) for correlation between removing the singularity in term C and improved convergence).
When the condensation rate based on SGR-P[0, 0, 0] is implemented in the CAM4 physics suite, most features of the simulated long-term climate remain similar to those in the original model based on M. Zhang et al. (2003) and Rasch and Kristjánsson (1998). These include the global-scale circulation and temperature pattern, as well as the geographical distribution and seasonal variations in rainfall and clouds (not shown). On the other hand, the cloud amount and the cloud radiative effects (CRE) show sizeable changes, as can be seen in Figure 3. The left column shows in black solid lines the 10-year annual mean, zonally averaged total cloud cover, long-wave (LW) and shortwave (SW) CRE in the two model configurations; the green and blue shading are the ±σ range where σ is the standard deviation of the 1-year averages. The right column shows the difference between each pair of 10-year averages in solid black lines and the ±σ range of the annual mean differences in light brown shading. The changes in these cloud-related quantities caused by switching to the SGR-P[0, 0, 0] based condensation rate are clearly above the natural variability. Specifically, larger cloud covers are seen at all latitudes and lead to an increase of the global mean from 0.54 to 0.77. The LW CRE also shows a systematic increase, with the strongest effects seen over the Southern Hemisphere storm tracks. The SW CRE is weakened in the low latitudes and strengthened in the middle and high latitudes. It is worth noting that such magnitudes of changes are not uncommon during model development, for example when a new parameterization is introduced or when uncertain model parameters are adjusted. The fact that changes in the short-term convergence behavior shown in Figure 2 can

Journal of Advances in Modeling Earth Systems
correlate to sizable changes in long-term climate suggests that the numerical aspects can be a substantial source of error and/or uncertainty in a global AGCM.
Although the convergence of the SGR-P[0, 0, 0] results is substantially better than the baseline model investigated by Wan et al. (2020) (see Figure 2), the results still show a degradation of convergence rate to 0.8 after 2 and 3 hr. This behavior is explored by ordering the grid cells according to the following measure: where T i ðt; ΔtÞ and T i ðt; 2ΔtÞ are the solutions for the ith grid cell using time step sizes Δt and 2Δt, respectively, and T * i ðtÞ is the reference solution defined in section 6.1. This measure is interpreted as how much Þ ), the contribution from the C term (ΔtC), and the changes in cloud liquid concentration caused by the resolved dynamics (ΔtA l ). The time series show that ΔtQ i ðtÞ is negative for most of the 3 hr simulation, even when q l i ðtÞ ¼ 0. This scenario activates a clipping limiter that prevents the condensation rate from driving the cloud liquid negative. The subgrid reconstruction profiles can be used to understand why the model is behaving this way.
After approximately 1,850 s, the mean condensation rate Q i ðtÞ in the problematic grid cell of Figure 4 becomes negative, meaning evaporation is predicted, primarily because of a negative C term and a zero B term (as A l is 0). That the cloud fraction is, however, monotonically increasing after 1,850 s is both counterintuitive and nonphysical. Recall that the C term for increasing cloud fraction (29) is derived by integrating the unaveraged condensation rate (14) over locations of cloud formation (28). As derived in section 3.2 and summarized in Table 1, cloud formation occurs at a location x over a time interval (t, t + τ) only when D(x, t, τ) ≥ D * (x, t), indicating there is sufficient moistening and/or cooling to overcome the saturation deficit. When SGR-P[0, 0, 0] is used, a nonnegative D( l n being greater than or equal to 1 Δt q sat n − q v n 1 − f ðt n Þ , which in turn translates to a nonnegative C term when A * l n ¼ 0 (as is the case in the problematic cell after 1,850 s.) A nonphysical, negative C term can result either from (i) incorrect assignment of a location to the cloud formation scenario (that does not meet the D ≥ D * criterion), assuming the atmospheric state at x is correct, or (ii) an inaccurate approximation of the atmospheric state at x that fails to meet the D ≥ D * criterion, assuming the determination of cloud formation scenario is correct. As mentioned in section 4.2, an empirical, RH-based expression of f is used to partition each grid cell in this work which means perspective (ii) is pertinent here. In summary, the nonphysical, negative C term in Figure 4 is a result from the piecewise-constantreconstruction SGR-P[0, 0, 0] not fulfilling D(x, t, τ) ≥ D * (x, t) in the cloud formation region predicted by the cloud fraction parameterization (21). Therefore the reconstruction profiles must be adjusted to avoid both poor convergence and nonphysical behavior.

Piecewise-Polynomial Reconstruction Results (SGR-P[1, 1, 1] & SGR-P[1, 1, 4])
The particular adjustment of the spatial profiles considered here to more easily meet the D ≥ D * criterion is to decrease the value of D * when f(t n ) < f(t n + 1 ) < 1/2, as is the case when the clipping limiter is active in the problematic grid cell of Figure 4. Recalling the one-dimensional grid cell x ∈ [0, L], the use of a piecewise-linear profile in place of the piecewise-constant profile for water vapor concentration will modify D * (x, t) in the cloud formation region (c(t n ) < x < c(t n + 1 )) by a factor of 2 Δt x − cðtÞ L − cðtÞ . Because f(t n + 1 ) < 1/2, one has that 2 x − cðtÞ L − cðtÞ < 1 for c(t n ) < x < c(t n + 1 ) and thus the factor reduces D * , making D ≥ D * more likely in the cloud formation region. With this motivation to use piecewise-linear profiles, the results for the piecewise-linear subgrid construction condensation rate (SGR-P[1, 1, 1]) are shown in Figure 5. The results show improved convergence through 3 hr, with a measured convergence rate of 1.0. Additionally, the problematic grid cell from Figure 4 no longer exhibits a strongly negative condensation rate that would drive the cloud liquid condensation negative (i.e., no strong clipping of Q i ðtÞ for t > 1850 s in ( Figure 5, right panel) As the simulation is carried out to 6, 9, and 12 hr, the convergence rate in Figure 5 drops to 0.9. To extend the improvement in convergence rate past 6 hr, an observation from Wan et al. (2020) is utilized. For the piecewise-constant subgrid reconstruction of A l (x, t) in SGR-P[0, 0, 0], the effect of the atmospheric flow (dynamical core) is distributed uniformly throughout the cell, as shown in the rightmost panel of Figure 1. When A l (x, t) < 0, this distribution implies that the atmospheric flow is drawing cloud liquid out of the grid cell equally from both the cloudy and cloud-free regions. It seems more physically valid for cloud liquid to be drawn out of the cloudy regions only; however, setting A l ðx; tÞ ¼ 0 in the cloud-free region of the grid cell breaks down when A l > 0 with f ¼ 0. While one could set A l ðx; tÞ ¼ 0 in the cloud-free region only when A l < 0, which is equivalent to what is done in Wan et al. (2020), it is not straightforward how to make the transition between the A l < 0 profile and the A l ≥ 0 profile smooth enough to be amenable to higher-order numerical techniques.

Journal of Advances in Modeling Earth Systems
With the piecewise linear subgrid reconstruction of A l (x, t) in SGR-P[1, 1, 1], Figure 1 illustrates how the effect of atmospheric flow is concentrated toward the cloudy region when that region exists. If p Al is increased past 1, the effect of the atmospheric flow is concentrated even more toward the cloudy region. This further concentration toward the cloudy region must be balanced with the smaller time steps required for the sharper spatial gradients. As such, a quartic profile is chosen for A l (x, t) resulting in a piecewise-linear-quartic subgrid reconstruction SGR-P[1, 1, 4]. Figure 6 shows that the linear-quartic profiles help extend first-order convergence further into the simulation, with measured convergence rate of 1.0 through 9 hr.
The climatological impacts of higher-order reconstruction are evaluated by conducting 10-year simulations similar to those presented in the previous subsection, using the CAM4 physics suite. In Figure 7, the zonally averaged annual mean cloud cover and CRE are compared between simulations using SGR-P[1, 1, 4] and SGR-P[0, 0, 0]. Again, changes that are sizable and significant compared with the uncertainties associated with natural variability are seen. There appear to be systematic decreases of cloud cover at all latitudes, with the largest changes occurring in the low latitudes and close to the Polar Regions. The changes in zonal mean CRE can be as large as about 10 W m 2 . Note that these changes appear to be smaller than the differences shown in Figure 3 between SGR-P[0, 0, 0] and the baseline formulation used in CAM4. This is consistent with the convergence testing results. It also suggests that the impacts of higher-order construction are smaller (albeit significant) than the impacts of changes in the condensation rate formulation resulting from the use of the SGR approach. Finally, it is worth noting that the SGR-P[1, 1, 4] based condensation rate does not significantly increase computational effort, costing 113.06 pe-hrs/simulated-year to produce 61.14 simulated-years/day compared to the baseline model cost of 113.87 pe-hrs/simulated-year to produce 60.70 simulated-years/day.

Conclusions
Direct simulation of the complex condensation and precipitation processes of water requires time scales and length scales beyond the capabilities of current computational resources available for global weather and climate models. Instead, modern global atmospheric models use parameterizations of the process known as large-scale condensation to describe the evolution of ensembles of cloud droplets and to take into account variations in the atmospheric state at spatial scales that are not resolvable by a global model. A parameterization for large-scale condensation can be derived from grid cell averaged values predicted by a global AGCM, as is the case with the parameterization studied by Wan et al. (2020). This approach is vulnerable to oversimplification that can lead to inconsistencies between the simulated cloud fraction and grid cell average cloud liquid concentration, which in turn lead to nonphysical values of in-cloud liquid concentration and cause time step convergence problems discussed in Wan et al. (2020). To avoid such nonphysical behaviors, an alternate parameterization is derived using a subgrid reconstruction methodology that respects subgrid physical relationships and constraints.

Summary
An unaveraged condensation rate was derived to describe the changes in atmospheric temperature, humidity, and cloud liquid concentration resulting from the phase changes of water. Scenarios of time evolution at a location are categorized according to whether the location is cloudy or cloud-free at the beginning and the end of an arbitrary time period. The criteria for this categorization, as well as the condensation rate in each scenario (where a negative value indicates evaporation), were derived by applying physically based constraints on the state variables. The connections between the unaveraged condensation rate in the unaveraged equations and the grid cell average equations in an AGCM were established using subgrid reconstruction and grid cell averaging.
The resulting formulation for the grid cell average condensation rate turned out to be similar to the one studied in Wan et al. (2020) but had important differences in the cloud growth/decay term. Specifically, the parameterization studied by Wan et al. (2020) (and the parameterization in CAM4 from M. Zhang et al. (2003)) followed Rasch and Kristjánsson (1998) in the use of an approximation that assumed (i) condensation in a cloud-growth (formation) scenario would work to bring the in-cloud liquid concentration to a value matching the existing clouds in the same grid cell, and (ii) evaporation in a cloud-decay (annihilation) scenario would work to remove the existing liquid water in that region. The new formulation derived here, in contrast, reveals that (i) externally induced moistening and/or cooling effects that exceed the saturation deficit in the cloud-growth (formation) region are offset by condensation, and (ii) the evaporation occurring in the cloud-decay (annihilation) region needs to remove not only the existing cloud liquid in the region but also eliminate any liquid newly imported by other processes (e.g., advection). Thus, the new formulation for the grid cell average condensation rate is a more complete and accurate description of the underlying physical processes. Additionally, it also inherently avoids the pathological situation of inconsistency that caused convergence problems in Wan et al. (2020), cf. formulation and discussions in section 4.2.
Evaluation of the different grid cell average condensation rates, which result from variants of the spatial reconstruction profiles used to connect unaveraged and grid cell average quantities, was based on the criteria of self-convergence properties and long-term climate impact. The first condensation rate evaluated is the one from choosing subgrid profiles that were piecewise-constant inside the cloudy or cloud-free areas for all the state variables and tendencies needed by the parameterization (SGR-P[0, 0, 0]). It showed both substantially improved convergence over the original parameterization studied in Wan et al. (2020), with first-order convergence achieved in 1 hr simulations, and significant long-term climate impact. The slightly degraded convergence rates in 2 and 3 hr simulations were found to result from the counterintuitive evaporation (instead of condensation) in cloud formation regions in some grid cells. Further analysis of both the numerical results and the derivation of the parameterization revealed an inconsistency in the reconstruct and average process.
For the purpose of avoiding highly complex model formulation and computationally expensive solution procedure, relative-humidity-based cloud fraction is used in place of a subgrid profile-based cloud fraction value to reconstruct and average the subgrid profiles. When applied to the piecewise-constant profiles of SGR-P[0, 0, 0], the averaging process used the relative-humidity based scheme to assign cloud formation to a region in violation of the scenario-categorization criterion. This led to the nonphysical evaporation that triggered strong clipping (limiter) to avoid negative liquid concentration, resulting in a degradation in convergence rate. Higher-order profiles that provide more subgrid variation were found to be useful in alleviating this consistency problem. For example, the piecewise-linear reconstruction (SGR-P[1, 1, 1]) avoids assigning cloud formation in violation of the scenario-categorization criterion and the corresponding strong clipper limiter activation, extending the first-order convergence to 3 hr. First-order convergence is extended even further by choosing a quartic profile for the cloud liquid tendency, giving the more physical effect of focusing the tendency in the cloudy region of the grid cell. The corresponding piecewise-linear-quarter reconstruction (SGR-P[1, 1, 4]) extends first-order convergence through 9 hr simulations while still having a significant impact on long-term climate.

Discussion
While this work shows the merits of a subgrid reconstruction methodology applied to only one parameterization, it serves as a demonstration of how a mathematically rigorous derivation can avoid nonphysical system states that would otherwise result from oversimplification and cause numerical problems. Clipping limiters can be implemented to abruptly force a nonphysical state back to a physically valid one. That said, such limiters can significantly affect convergence properties of the numerical solutions. Furthermore, the use of such limiters leaves the root causes of the nonphysical behavior unaddressed, which can have substantial impact on the long-term climate. If instead, the physical constraints are respected and maintained throughout the parameterization derivation, such as cloud liquid presence requiring vapor saturation, no clipping limiter or other abrupt adjustment should be needed, because the numerical solutions will behave in a correct way for a correct reason.
This work also provides a deep understanding of the cause of pathological behavior resulting from inconsistencies between the choices of cloud fraction scheme and reconstruction profiles. This understanding has implications beyond the relatively simple parameterization investigated here. In atmospheric models, it is common to have cases where separate equations are formulated-and numerically solved-for closely related quantities. Physically meaningful and important relationships between such quantities can get lost during the simplification of model equations or as a consequence of discretization. One such example is the advection of linearly correlated trace species. The loss of correlation due to discretization caught a substantial amount of attention in the past two decades, first among developers of chemical transport models and later in the wider atmospheric modeling community (see, e.g., Jöckel et al., 2001;Lin & Rood, 1996;Thuburn et al., 2014;Zhang et al., 2008). By now, preservation of linear correlation has become one of the basic requirements for modern advection schemes (see, e.g., Bosler et al., 2019;Lauritzen et al., 2011;Ullrich et al., 2013).
Another example of a frequently encountered inconsistency can be found between particle mass and number concentrations, such as those of cloud droplets or aerosol population. Because the mass and number concentrations can be affected by different physical processes whose representation in an AGCM can involve different levels of simplification, and because the mass and number concentrations are often transported as separate tracers in modern model, it is not unusual to see nonphysical situations where there is zero mass concentration with nonzero number concentration, or vise versa. In E3SM and its recent predecessors, for example, such cases are handled by setting both concentrations to 0. The impacts of the inconsistencies behind such clipping are not yet clear.
A third example the authors have run into is with the prediction of the second and third moments of the subgrid probability distribution of vertical velocity in the turbulence and clouds parameterization CLUBB (Cloud Layers Unified By Binormals, Golaz et al., 2002;Larson & Golaz, 2005). Because these moments are calculated by solving separate partial differential equations, one can find situations where there is zero second moment with nonzero third moment. Note that the relative magnitudes of the two moments determine the skewness of the subgrid statistical distribution of the vertical velocity, which then has a direct impact on the transition between stratocumulus clouds and shallow convection. This transition is a well-known challenge in parameterization development. The impact of inconsistencies related to model formulation or numerical methods on the cloud properties simulated by CLUBB, and by the global models in which CLUBB has been embedded (e.g., Bogenschutz et al., 2013;Guo et al., 2014;Rasch et al., 2019), is also unclear.
From the results presented in this paper, the authors speculate that investigations into the above mentioned inconsistency issues will lead to opportunities to further improve the related parameterizations. Modern parameterizations of turbulence, convection and cloud/aerosol microphysics are very complex; simplifications of the equations and approximations in the numerical methods are necessary to keep the computation tractable. While physical intuition is important in guiding such simplifications and approximations, mathematically rigorous derivation can also provide useful insights. In addition, as seen in this work, the use of flexible mathematical frameworks, including the subgrid reconstruction methodology used here, can provide opportunities to improve the inherent consistency of the simulations.
where δ is the one-sided Dirac delta function and t CA 1 ðxÞ, t CA 2 ðxÞ, … are all the times that the cloud annihilation scenario occurs at location x. Similar steps will give the following for the cloud formation scenario at location x and time t: where t CN 1 ðxÞ, t CN 2 ðxÞ, … are all the times that the cloud formation scenario occurs at location x.

Data Availability Statement
The code used to produce the long-term climate results is provided in Vogl et al. (2020a). The code used to produce the self-convergence test data and figures is provided in Vogl et al. (2020b).