## Abstract

In this paper a large-eddy “smooth” cloud (SC) model will be presented with smooth implying that the entire model converges under a Newton-based solution procedure or that time scales within the SC model are being resolved. Besides ensuring that time scales within microphysical parameterizations are resolved, convergence of Newton’s method requires that advection schemes near cloud boundaries should not induce fast time scales. For example, flux-corrected transport (FCT) schemes that force cloud variables to stay oscillation free near boundaries are typically not differentiable in time and hence may prevent convergence of Newton’s method. To circumvent the use of a FCT scheme, an alternative approach, a cloud-edge (CE) diffusion-based approach, will be presented in this paper. Since the diffusion produced by the CE approach could conceivably lead to the fictitious evaporation of a real cloud, the first major point of this paper will be to document that the SC model when employing an evaporative limiter is able, like most traditional large-eddy cloud models, to reasonably reproduce nondrizzling stratus clouds observed during flight 1 of the Second Dynamics and Chemistry of Marine Stratocumulus field study (DYCOMS-II). However, the SC model obtains the accuracy offered by higher-order time-stepping approaches, unlike most traditional cloud models. In fact, temporal errors from the SC model are shown to be at least two orders of magnitude smaller than those of a traditional large-eddy cloud model. Hence, the second major point of this paper will be to demonstrate the consequence of these large temporal errors found in traditional large-eddy cloud models, that is, the inability to accurately track an identifiable cloud feature in time.

## 1. Introduction

In Reisner et al. (2005) a numerical framework based upon the Jacobian-free Newton–Krylov method (JFNK; see Table 1 for a list of key abbreviations used in the paper) was presented that enabled second-order-in-time accuracy to be achieved for both a highly idealized moist bubble and a hurricane. Though the JFNK method was shown to produce much smaller time errors than a leapfrog semi-implicit (SI) numerical approach, the JFNK method unfortunately employed a large constant eddy diffusivity to resolve cloud edges (CEs) that could, in certain situations, lead to the complete evaporation of certain types of clouds (e.g., stratus clouds). This paper presents an alternative applied first within large-eddy simulations of stratus clouds.

Though the proposed alternative is conceptually very simple with high levels of diffusion being applied only near cloud edges within a smooth cloud (SC) model, the impact of this differential diffusion with regard to a model’s ability to accurately represent a realistic cloud field, such as a stratus deck, is not clear. For instance, even though the diffusion is only active near the edges, the diffusion could, in conjunction with evaporation, still lead to a fictitious thinning of a cloud feature. To help prevent this thinning, an evaporative limiter was implemented that minimizes evaporation when the time scale of the extra diffusion exceeds the time scale of diffusion associated with a traditional turbulence model.

Because the additional diffusion produced by the CE approach could lead to the spurious modification of cloud features, the first major point of this paper will be to demonstrate the ability of the SC model to reproduce a realistic cloud field as well as a traditional large-eddy cloud (TC) model. To address the first point both the SC model and a TC model will be used in large-eddy simulations of stratus clouds observed during the Second Dynamics and Chemistry of Marine Stratocumulus field study (DYCOMS-II; Stevens et al. 2005). Note, in this paper a TC modeling approach is defined to be a numerical procedure employing two time-split solution approaches, a process to advance the cloud model in time and step over fast time scales, such as the SI procedure, and an advective operation, such as flux-corrected transport (FCT; Zalesak 1979) that selectively adds sufficient diffusion in time to minimize numerical oscillations near cloud edges. Additionally, because a majority of large-eddy cloud models used within the atmospheric community to simulate stratus clouds use one or both time-split procedures (i.e., see appendix A of Stevens et al. 2005), it is hypothesized that any results from the TC modeling approach used in this paper will be similar to results produced by any one of those models.

Since most TC models employ FCT procedures (Skamarock 2006) to implicitly resolve cloud edges and being that these numerical procedures either are not differentiable or at the very least induce very fast time scales into the discrete equation set, the second major point of this paper will be to illustrate the negative consequence—in terms of the production of large temporal errors—of employing nearly discontinuous numerical procedures within TC models. In contrast, since the discrete SC model is able to achieve the accuracy offered by higher-order time-stepping procedures, temporal errors within the SC model framework are shown to be much smaller than the errors produced by TC models. Furthermore, numerical problems associated with not resolving time and spatial scales near cloud edges have been recognized since the onset of cloud modeling with numerous papers (Klaassen and Clark 1985; Grabowski and Smolarkiewicz 1990; Stevens et al. 1996; Margolin et al. 1997; Jeffery and Reisner 2006; Grabowski and Morrison 2008; Morrison and Grabowski 2008) either discussing the problems or presenting possible solutions to the problems. The SC modeling approach represents another possible solution to minimize numerical issues near cloud boundaries.

Two particular items distinguish the current SC modeling approach from other traditional approaches. First, it is required that the entire equation set converge within the JFNK solution procedure, thus mandating within the discrete model that the smallest values for cloud quantities are near the Newton tolerance. But, the small background values have little impact outside of cloud boundaries with these values being subtracted from forcing terms such as buoyancy, and therefore the small values simple serve to ensure convergence of the discrete model to the specified nonlinear tolerance, *ε*_{nl}. Likewise, since advection in the JFNK cloud model is implicit with both the advected quantity and the velocity being at the new time level, the accuracy of this nonlinear term and the models ability to keep cloud variables positive is directly a function of the magnitude of the nonlinear tolerance.

Keeping in mind the inherent limitations of employing a continuous approach near cloud interfaces instead of a discrete particle-based approach (Andrejczuk et al. 2008), which can resolve the subgrid movement of cloud boundaries, the second major item that distinguishes the current approach from other approaches is that cloud edges are resolved over more than one grid cell. By explicitly adding diffusion and enabling cloud fields to be resolved near boundaries, oscillations associated with higher-order advection schemes can be minimized and high temporal accuracy concerning the movement of clouds can be achieved. Furthermore, because the FCT procedure typically reverts to the highly diffusive upwind advection scheme to resolve cloud boundaries, this procedure can *sporadically* add as much or even more diffusion than the CE approach. And, since the FCT procedure cannot resolve the subgrid movement of cloud boundaries (Margolin et al. 1997), oscillations are produced in both temperature and water vapor fields that can than feed back into cloud fields. Though these oscillations can be somewhat minimized by using a prognostic supersaturation equation (Grabowski and Morrison 2008) instead of the water vapor equation, the inherent numerical problems associated with employing a FCT-based advection procedure still remain.

The obvious trade-off in the current approach with regard to traditional approaches is that more than one grid cell is required to resolve a cloud edge and hence this realization implies that the cloud feature of interest should be present over at least a few grid cells. Likewise, the other major limitation of the approach is the introduction of a tuning coefficient within the evaporative limiter. But, since little control exists with regard to the amount of CE evaporation produced by TC models, at least the SC model allows for the possibility to limit excessive evaporation near cloud edges.

The remainder of the paper demonstrates the two major points of this paper with the next section presenting the analytical equations of the cloud model; section 3 presents simulations of the stratus deck observed during DYCOMS-II used in validating the SC model; section 4 illustrates the consequences of large temporal errors in two-dimensional (2D) simulations of a moist bubble interacting with the stratus deck observed during DYCOMS-II; and the final section offers a brief summary and future directions.

## 2. Cloud model

In this section the analytical equations for the cloud model are presented. The dynamical equation set involves expressions for the conservation of momentum, energy, and mass of the gas, whereas the bulk microphysical model that is appropriate for a nonprecipitating stratus cloud is composed of conservation equations for cloud number concentration and cloud water density. Smooth expressions that allow for momentum, energy, and mass exchanges to occur between the dynamical and microphysical models are also included in the cloud model (presented in the appendix as well as a turbulent kinetic energy (TKE) model to simulate subgrid fluctuations in dynamical and cloud fields.

### a. Analytical equation set

Even though the discrete model is formulated in a three-dimensional (3D) generalized coordinate frame, for ease of presentation the equation set will be represented in 3D Cartesian space.

In this framework the momentum equations are expressed as follows

where all indices range from 1 to 3; *u*^{1′} and *u*^{2′} are the Cartesian gas velocities in the horizontal, *x*^{1′} and *x*^{2′}, directions; *u*^{3′} is the Cartesian gas velocity in the vertical, *x*^{3′}, direction; Ω_{1′} = 0 with Ω_{2′} = 2Ω sin*ϕ* and Ω_{3′} = 2Ωcos*ϕ* being the *x*^{3′} and *x*^{2′} components of the earth’s rotation axis at the latitude *ϕ*; *ρ* is the gas density of the air, *ρ* = *ρ _{d}* +

*ρ*with

_{υ}*ρ*the dry air density and

_{d}*ρ*the water vapor density;

_{υ}*p*′ =

*p*−

*p*is the pressure perturbation with

_{e}*p*the pressure of the gas computed using the ideal gas relationship [i.e., Eq. (11) in Reisner et al. 2005)] and

*p*=

_{e}*p*(

_{e}*z*) the environmental pressure;

*ρ*′ =

*ρ*−

*ρ*is the density perturbation with

_{e}*ρ*=

_{e}*ρ*(

_{e}*z*) the environmental density and the buoyancy force only being applied in the vertical direction,

*δ*

_{i′3};

*ρ*′ =

_{c}*q*−

_{c}ρ*ρ*is the perturbation cloud water density with

_{ce}*q*as the “specific” cloud water mixing ratio [i.e., the mass of cloud water per the mass of the gas (water vapor plus dry air)] and

_{c}*ρ*=

_{ce}*q*, where

_{ce}ρ*q*= 0.1

_{ce}*ε*

_{nl }

*q*is the environmental cloud water mixing ratio with

_{cb}*q*= 1 kg kg

_{cb}^{−1}a constant base-state

*q*field;

_{c}*u*

_{e}^{i′}are the environmental winds in all three directions;

*τ*

^{i′j′}= (∂

*u*

^{i′}/∂

*x*

^{j′}) + (∂

*u*

^{j′}/∂

*x*

^{i′}) − (2/3)

*δ*

^{i′j′}∂

*u*

^{s′}/∂

*x*

^{s′}is the strain-rate tensor; and the parameters

*δ*and

*ε*are used to simplify the writing of the gravitational and Coriolis terms, respectively (see chapter 2 of Pielke 1984). Note that constants such as the earth’s gravity

*g*are defined in Table 2 and key variables within the cloud model such as

*u*

^{i′}are defined in Table 3.

The energy equation is expressed as

where *θ* is the potential temperature *θ* = *T*(*p _{o}*/

*p*)

^{Rd/Cp}, with

*T*the temperature of the gas,

*f*

_{energy}= (

*θρL*/

*TC*)

_{p}*f*

_{mass}represents the release of energy during phase conversion of a given amount of mass

*f*

_{mass}, associated with various microphysical source terms (described in the appendix);

*f*

_{rad}is a simple parameterization for cloud-top longwave-induced cooling (Stevens et al. 2005); and the diffusional flux of potential temperature fluctuations being defined in the

*x*

^{1′}direction is

*F*

_{θ}^{1′}=

*ρκ*(∂

*θ*/∂

*x*

^{1′}).

The diffusion coefficient *κ*, found in Eqs. (1) and (2), as well as all other scalar equations, is determined from a TKE model with this equation being expressed as

with the first term on the right-hand side being the shear generation of turbulence, the second term the buoyancy generation of turbulence, the third term the dissipation of turbulence, and the last term the diffusion of turbulence where *κ* = 0.09*L _{s}*TKE with

*L*being the turbulence length scale.

_{s}In contrast to the momentum and energy equations, the diffusion coefficient found within the turbulent fluxes for all other predictive equations is represented as follows:

where with from the CE approach, and *ψ* represents the quantity being diffused [e.g., *ψ* = TKE in the turbulent fluxes of Eq. (3)].

The conservation equation for water vapor density *ρ _{υ}* can be written as follows:

where *ρ _{υ}* =

*ρq*with

_{υ}*q*the specific humidity. Note, the conservation equation for

_{υ}*ρ*is nearly identical to Eq. (5), but that the turbulence term is not present. Likewise, the conservation equation for

*ρ*is similar to Eq. (5), aside from the term

_{c}*f*

_{mass}being positive and the addition of a term to represent the fall of

*ρ*.

_{c}The conservation equation for cloud number density *N _{c}* can be written as follows:

where the fall term fall* _{Nc}*and the source term

*f*

_{number}are described in the appendix.

### b. CE approach

Using dimensional analysis to assist in the formulation of the diffusivity coefficient found in the CE approach, the extra diffusion for *ρ _{c}* (all other scalar fields are treated in a similar manner) is expressed as follows:

where is a closure coefficient; the absolute value of a velocity component |*u*^{i′}| is approximated during a simulation by the following smooth function, |*u*^{i′}| ≈ *u*^{i′} tanh(100*u*^{i′}); *L*_{extra} is a length scale of the extra diffusion that is set to Δ*x*^{i′} during a simulation; *ε* = 10^{−12} kg kg^{−1} is a constant to prevent division by zero; and Δ_{i′}*q _{c}* represents the change of the

*q*field in the

_{c}*x*

^{i′}direction. To ensure that the cloud-edge diffusion coefficient does not introduce fast time scales, the coefficient has been limited as follows:

where

with Δ*t*_{ref} = 10Δ*t*_{sound} and Δ*t*_{sound} being the sound wave time scale.

In this paper, two approaches for determining will be utilized, with the first approach employing a constant and the second approach relating to an advective operator that determines whether or not *ρ _{c}* remains positive. For the first approach, simple scaling arguments along with some trial and error were used to deduce that ≈ 0.01. Obviously, could be increased over this value; however, in this context becomes more of a traditional turbulence tuning parameter and less of a means to ensure that the higher-order interpolation employed in the advection scheme does not produce negative cloud values near cloud boundaries.

The second approach is a two-step process with the first step in discrete notation being as follows:

and the second step is

where *ϕ*^{CE−ADV} is a parameter controlling how quickly transitions between 0 and 1; * indicates any time level between the old and new time level or a pseudo time level; *ũ*^{i′} represents normalized cell-face velocities (described in section 2c) with *ρ̃* being interpolated to a given cell face by an advective scheme, and the implied spatial subscript *i* denotes the locations *x ^{i}* =

*i*Δ

*x*on the grid. Notice, as

^{i}*ϕ*

^{CE−ADV}increases in magnitude, tanh[

*ϕ*

^{CE−ADV}(

*ρ*

_{c}

^{1}/

*ρ**)] only occasionally becomes negative at cloud edges, and hence the application of diffusion becomes highly sporadic in time and essentially mimics what is being done by the FCT procedure at the cloud edge (i.e., occasionally adding diffusion via the upwind advection scheme). Additionally, since the amount of diffusion being added is a known quantity, any evaporation brought on by this diffusion has been limited by an evaporative limiter with details concerning this limiter presented in the appendix.

_{c}### c. Discrete model

A variety of time-stepping approaches have been utilized within the SC model and can be symbolically represented as follows:

where **x** is a vector containing all variables of interest; *n* + 1 or *n* represents either the new or old time levels; 𝒫 denotes some time-stepping procedure or process [i.e., the SI process, see Reisner et al. (2005)]; and **F** is a vector representing all the equations contained within a cloud model. Thus, in practice, all time-stepping procedures utilize the same equation set and hence the impact of these different time-stepping procedures can be readily investigated.

Excluding the momentum equations, the discrete form of all other predictive equations can be represented as follows:

where a coincident or nonstaggered grid was employed and ℳ denotes terms associated with the microphysical model (see the appendix). Additionally, in Eq. (13) all quantities not at the new or old time level are implied to be at time level *. In Eq. (13), 𝒫 represents a variety of time-stepping procedures with a key aspect to ensure efficiency of the implicit procedures being the use of a physics-based preconditioner (Reisner et al. 2005), which is also used as the SI solution procedure. Furthermore, when the discrete equation set is iterated the chosen nonlinear convergence criteria is as follows:

where *ε*_{nl} = 1 × 10^{−5} and the norm in Eq. (14) is defined by with *n*_{tot} as the total number of relevant grid variables and *k̂* as the iteration index.

The advective quantities found in Eq. (13) are composed of a cell-face Courant number:

and a quantity that is interpolated to the appropriate face, that is, , via the quadratic upstream interpolation for convective kinematics advection scheme (QUICK; Leonard and Drummond 1995) or the QUICK scheme including estimated streaming terms (QUICKEST) with these quantities have the possibility of being limited by a FCT procedure (Zalesak 1979). The choice of two advection schemes, instead of one, was mandated by the requirement that the advection scheme for the SC model employing higher-order-in-time solution procedures be relatively smooth and produce temporal errors that are clearly better than first order, whereas for traditional approaches the primary requirements are the stability and efficiency of the procedure (i.e., the QUICK scheme is not stable when using a one-step time-advancement procedure). Also, averaging the cell-face Courant velocities in time as well as in space to form a pseudo “C grid” reduces numerical issues associated with a collocated mesh or “A grid” and helps improve model stability (Reisner et al. 2000; Smolarkiewicz and Margolin 1998).

The discrete form of the momentum equations is nearly identical to Eq. (13), but with the inclusion of the pressure gradient terms that employ a second-order-in-space central difference approximation (except at the top and bottom boundaries) and the Coriolis terms calculated at the appropriate time level. Also, the diffusion terms in the momentum equations are calculated in the same manner as in the scalar equations with each component of the stress tensor being computed on the appropriate cell face and then differenced across the grid cell.

## 3. DYCOMS-II simulations

Given that inherent observational uncertainties and the idealizations of any numerical model make quantitative model error assessment difficult to formulate, the primary goal of this section will be to demonstrate that the SC modeling framework produces results that are both comparable to the TC modeling framework employed in this paper and to the other TC approaches shown in (Stevens et al. 2005). Likewise, noting the delicate balance between longwave cooling, subsidence, surface fluxes, condensation, and advection that is required to maintain a stratus deck, it is not obvious whether the additional diffusion being added by the CE approach will upset this balance or the evaporative limiter will be robust enough to prevent too much thinning of the cloud deck. Additionally, one-dimensional tests employing both discrete frameworks have been conducted to document that they do indeed converge to approximately the same solution as the grid spacing goes to zero. However, given the rather cumbersome nature of one-dimensional tests (i.e., they mainly measure linear advective errors and not errors from nonlinear feedbacks), results from these tests have not been included. This section serves as the primary documentation that the SC model produces spatial errors that are comparable to those produced by the TC model.

### a. Model setup

Results from four simulations, SC1, SC2, SC3, and TC1, are presented (see Table 4 for an overview of the setup). Specifically, the setup for the four simulations are as follows: simulation 1, SC1, employed the SC model with = 0.01 for *q _{υ}*,

*q*,

_{c}*N*, and TKE ( = 0 for the other equations); simulation 2, SC2, was identical to SC1 except the evaporative limiter was active; simulation 3, SC3, employed the SC model and the second advective procedure—designed to mimic the FCT scheme—to determine with

_{c}*ϕ*

^{CE−ADV}= 1 for

*q*,

_{υ}*N*, and TKE and

_{c}*ϕ*

^{CE−ADV}= 100 for

*q*; simulation 4, TC1, utilized the FCT procedure instead of the CE approach to ensure cloud variables remain positive. Simulations SC1 and SC2 employed the QUICK advective scheme, whereas simulations SC3 and TC1 employed the smoother variant of the QUICK scheme, QUICKEST, for advection.

_{c}For efficiency, simulations SC1 and SC2 employed an iterative “Jacobi like” SI solution procedure to ensure convergence to the specified nonlinear tolerance; that is, the cost of inverting the entire Jacobian matrix [e.g., Eq. (21) in Reisner et al. (2005)] is considerably higher than inverting a single scalar pressure equation [e.g., Eq. (A39) in Reisner et al. (2005)] for each nonlinear iteration. Additionally, both simulations SC1 and SC2 utilized the second-order-in-time Crank–Nicholson (CN) time-stepping procedure. Because the numerical formulations of either simulation SC3 or TC1 were not differentiable in time and produce a relatively large nonlinear tolerance regardless of the number of SI iterations, the two simulations employed a standard one-step SI procedure with a first-order-in-time Euler time-stepping approach.

The initialization of model variables for the four simulations closely follows what was described in Stevens et al. (2005) with mean fields using data and figures contained within section 2a(1) of that paper, except for the *N _{c}* field that was initialized from data shown in Fig. 2 (

*N*≈ 100 cm

_{c}^{−3}) of Twohy et al. (2005). To initiate vertical motions, a weak moist bubble was introduced into the model domain of the following form:

where *f*_{time} = (*t* − 330 s)/30 s; and *f*_{space} = dis(*x*^{i′})/180 m with dis(*x*^{i′}) the distance in meters from the center of the bubble (*x*^{1′} = 1600 m, *x*^{2′} = 1600 m, and *x*^{3′} = 260 m). Note, the small magnitude of the source function was by design; a larger value results in an unacceptable vertical rise of the stratus deck.

A variable vertical resolution using 115 grid points produced by the following relationship (*x*^{3′} is in meters):

was used in a majority of the simulations along with a constant horizontal resolution of 40 m over 80 × 80 grid points. All simulations were run for a time period of 4 h and to compare against the averaged data presented in Stevens et al. (2005), key model quantities, such as *q _{c}*, were averaged in space and time during the last 1 h of model integration.

### b. Results

Figure 1 shows 3D isosurfaces from the four simulations with simulation SC1 the only simulation of the four that produces a somewhat broken-up cloud deck. Plots of the mean vertical profiles of *q _{c}* and

*N*(see Figs. 2a,b) illustrate that the lack of cloud water evident in Fig. 1 for simulation SC1 is also apparent in the mean profiles with simulation SC1 producing a considerably thinner cloud than the other simulations. Likewise, Figs. 2a,b clearly reveal the impact of the evaporative limiter with simulation SC2 producing the thickest stratus deck with the cloud bottom from simulation SC2 being considerably lower and closer to the observed cloud bottom (dashed black line in Fig. 2a) than what is produced by the other three simulations.

_{c}Given that simulations SC3 and TC1 were of similar design, it is not surprising that their mean *q _{c}* and

*N*fields (see Figs. 2a,b) are roughly similar; however, there are significant differences in mean vertical motion statistics (see Figs. 2c,d) from the two simulations with both the variance and the third moment of the vertical motion field being significantly larger in simulation SC3, but still less than the observed variance shown in Stevens et al. (2005, see their Fig. 5). This result suggests that evaporation-induced entrainment events are significantly stronger in simulation SC3 than simulation TC1 with analysis of both simulations showing that the diffusion produced by the CE approach in simulation SC3 is rarely active, but in simulation TC1 the FCT procedure occasionally is activated leading to diffusion and a reduction in the intensity of the simulated entrainment events. Likewise, Figs. 2c,d show that in simulations SC1 or SC2 the large amount of diffusion or the lack of evaporation found in the respective simulations also leads to a significant lessening in mean vertical motion statistics.

_{c}A spatial analysis (not shown) of the evaporative limiter within simulation SC2 reveals that this limiting primarily occurs at cloud top or where the highest gradients in *q _{c}* exist. Furthermore, the spatial analysis reveals that the majority of the limiting is associated with the relatively large amounts of cloud-edge diffusion being added in the horizontal direction, instead of the vertical, due to the presence of higher wind speeds in the horizontal versus vertical direction. Additionally, evaporative limiting not only leads to more

*q*at cloud top, but because of the falling term [e.g., Eq. (A20)], this greater amount of

_{c}*q*falls downward, thus producing a thicker cloud.

_{c}Given the rather large impact of the evaporative limiter on the dynamical fields, an additional sensitivity simulation was conducted in which the value of the tuning parameter, *ϕ*^{evtune}, found within the evaporative limiter was reduced from 10 to 1. Figures 3a–d shows that while this reduction has only a small impact on mean *q _{c}* and

*Nc*fields, the reduction does have a rather large effect on the vertical motion statistics with a significant increase in vertical velocity statistics, especially in the variance. Hence, this result tentatively suggests that the optimal value of

*ϕ*

^{evtune}should be between 1 and 10.

Overall, this section clearly showed that the SC modeling framework did not produce results outside the realm of observations, the TC model, or what was produced by other large-eddy TC models shown in Stevens et al. (2005). Furthermore, even though the FCT procedure appears to produce less diffusion than the CE approach employing a constant , the comparison of simulation SC3 against simulation TC1 did demonstrate that the FCT procedure does introduce considerable diffusion with the magnitude of this diffusion being reasonably close to what is produced by simulation SC3. Unfortunately, as will be shown in the next section, the sporadic application of diffusion in simulations SC3 and TC1 produce large temporal errors that prevent either numerical procedure from accurately tracking features in time.

## 4. Temporal error

As mentioned in the introduction, the primary point of this section is to document that the nearly discontinuous-in-time FCT procedure introduces very fast time scales into the discrete equation set that result in the production of large temporal errors, as opposed to the smooth CE approach that induces only small temporal errors. As in section 2 of Mousseau et al. (2002), temporal error could be investigated in an analytical setting by introducing Taylor series expansions of the following form:

into the discrete equations of the cloud model. Unfortunately, undertaking this procedure within the TC modeling framework is difficult, if not impossible, because of the sporadic on-and-off switching of the FCT procedure that is difficult to capture in an analytical setting. Likewise, it is not obvious whether this switching introduces a zeroth-order temporal error into the discrete equation set [i.e., temporal derivatives in Eq. (18) are nearly independent of time step size].

To demonstrate in a computational setting that large temporal errors are produced by the nearly discontinuous-in-time FCT procedure, simulations of a moist bubble impacting the stratus deck observed during DYCOMS-II (Stevens et al. 2005) will be conducted with the temporal error being computed by comparing results from simulations employing various time step sizes against reference simulations using very small time steps. Because of the computational expense of computing the reference simulations, simulations to be presented in this section will be two-dimensional. For an explanation as to why the numerical tests conducted in this section primarily address temporal error see Figs. 3–4 of Reisner et al. (2003) along with the corresponding discussions. Furthermore, though this section is intended to illustrate the negative consequences of employing the FCT procedure in place of the smooth CE approach, it should be stressed, once again, that any numerical procedure that produces either an ill-conditioned Jacobian matrix and/or prevents Newton’s method from converging should induce very large temporal errors.

### a. Model setup

As in the previous section individual simulations will be denoted as SC1, SC2, SC3, and TC1; however, unlike the previous section, simulations will be grouped into sets with the only difference between the various simulations within each set being the time step size (i.e., SC1–1/16 denotes a simulation within set SC1 that employed a 0.0625-s time step; see Table 5 for an overview of the current setup). Also, unlike in the previous section, a 2D version of the cloud model will be utilized with all simulations employing 160 horizontal grid points and 64 vertical grid points with a constant resolution of 20 m. But, initialization of model variables is the same as in the previous section.

To compute temporal error for simulation sets SC1 and SC2 that utilized a fourth-order explicit Runge–Kutta time-stepping procedure, reference simulations employing a time step size of 0.1Δ*t*_{sound} ≈ 0.006 25 = 1/160 s were run using a three-stage semi-implicit Runge–Kutta procedure (e.g., SC1–1/160r with “r” denoting a reference simulation for SC1), whereas temporal error calculations for simulation sets SC3 and TC1 were calculated by comparing against reference simulations using a time step size of 1/1600 s with these simulations employing the SI solution procedure along with the first-order-in-time Euler time-stepping approach. For each reference simulation, a smooth water vapor source function of the following form:

was used to initialize a moist bubble where *f*_{time} = (*t* − 330 s)/30 s; and *f*_{space} = dis(*x*^{1′}, *x*^{3′})/180 m with dis(*x*^{1′}, *x*^{3′}) the distance in meters from the center of the bubble (*x*^{1′} = 900 m, *x*^{3′} = 260 m). All simulations were then initialized from output obtained from the reference simulations at *t* = 360 s with these simulations continuing for another 1200 s. The use of either the explicit Runge–Kutta approach in simulation sets SC1 and SC2 or the SI procedure in simulation sets SC3 and TC1 was based on numerous sensitivity tests examining temporal accuracy and efficiency of various time-stepping procedures (i.e., CN time stepping within the JFNK procedure), with the two chosen methods minimizing the product of the two quantities. Furthermore, as noted in the previous section, it makes little sense to employ higher-order time-stepping procedures within the TC model when a numerical method, such as FCT, is at best first-order-in-time accurate.

The error measure (e.g., in the *q _{c}* field) used to quantify differences between a given simulation and a reference simulation is defined as follows:

with q_{c}^{ref} coming from a reference simulation and *q _{c}^{l}* being the computed value at a grid point

*l*.

### b. Results

Figure 4 shows a time history of the moist bubble, as indicated in the *q _{c}* field, from the reference simulation, SC2-l/160r, with this figure demonstrating the impact of the bubble with the stratus field along with the development of very sharp cloud boundaries around the bubble’s top and sides. Figures 5 –6 show

*q*and wind vector fields along with differences in these fields with regard to their respective reference simulations at the end of simulations SC2–1/8 and SC2–1/160, and simulations TC1–1/8 and TC1–1/160. These two figures clearly reveal the inability of the traditional approach to accurately simulate the interaction of the moist bubble with the stratus deck, whereas, in stark contrast, the SC model is able to accurately simulate the temporal evolution with regard to the reference simulation employing a very small time step.

_{c}The large visual differences between the two approaches also translates into large differences in temporal error with Fig. 7 showing at least a three order difference in temporal error between simulations within set SC2 versus set TC1. Likewise, Fig. 7a clearly shows that all simulations within set TC1 experience a rapid error growth within the first 60 s with the errors being roughly the same regardless of time step size. But, in marked disagreement with results from simulation set TC1, the errors (see Fig. 7b) from simulations within set SC2 do show a rapid decrease in error with decreasing time step size and while the error growth does increase during the simulations, it is not as rapid as simulation set TC1 or set SC3 (not shown).

Table 6 further quantifies the errors from the four simulation sets as a function of time step and illustrates that the errors from simulation sets TC1 and SC3 are many orders of magnitude larger than from either simulation set SC1 or SC2. Once again, the lack of temporal error convergence is clearly evident in the table with simulation set SC3 producing errors that are even larger than simulation set TC1 with the fast on–off switching of the CE diffusion found within the second advective approach of the CE approach being responsible for these errors. Hence, Table 6 and the figures shown in this section together illustrate the negative consequences of using a numerical procedure that is nearly discontinuous in time; the generation of large temporal errors that do not substantially decrease with time step size and prevent a key feature, i.e., the moist bubble, from being tracked accurately in time.

## 5. Conclusions and future directions

In this paper an SC model was presented that enables nonlinear convergence in time for problems containing relatively sharp gradients in cloud variables. A key aspect to achieve nonlinear convergence was the CE approach that replaces the nearly discontinuous in time FCT advection procedure found within TC models. To demonstrate in a realistic setting the fidelity of the CE approach with respect to the FCT procedure, simulations of stratus clouds observed during DYCOMS-II (Stevens et al. 2005) were conducted with the SC model producing results comparable to both observations and what was produced by a TC model. But, as shown in the previous section, the chief benefit of the SC model is its ability to produce very small temporal errors that rapidly decrease with time step size.

The smaller temporal errors of the SC model were the result of employing a discrete representation of the analytical equation set that was differentiable and could realize the accuracy offered by higher-order temporal differencing schemes (i.e., almost all flux-limited advection procedures that guarantee monotonicity are not differentiable in time). But, though the CE approach has shown the potential to minimize numerical errors found near cloud boundaries, neither the approach of applying differential diffusion nor the evaporative limiter is unique with future work still being required to refine these procedures. For example, future questions that should be addressed are whether or not *ϕ*^{evtune} varies with each condensate species and what is the optimal value for *ϕ*^{evtune}. Another point not shown in this paper is the impact of the CE approach on spatial errors with preliminary analysis of much higher-resolution moist bubble simulations suggesting that the SC model, because of its use of the evaporative limiter within the CE approach, is able to produce results at coarser resolutions that are comparable to higher-resolution findings.

A future avenue would be to embed all or part (just the CE approach) of the SC modeling approach within an operational weather prediction model to help determine whether or not the approach being advocated in this paper improves model predictability. While the CE approach would be relatively easy to incorporate into an existing model framework, the most difficult part would be to ensure that time scales in physical models are properly limited such that they do not induce rapid dynamical changes in predicted quantities. Otherwise, the time step size may need to be made unreasonably small before an examination on the impact of the CE approach on model predictability can be made.

## Acknowledgments

This research was supported by Los Alamos National Laboratory’s Directed Research and Development Project “Resolving the Aerosol-Climate-Water Puzzle (20050014DR)” with Dubey Manvendra serving as principal investigator for this project. The Los Alamos National Laboratory was sponsored under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy under DOE Contract W-7405-ENG-36. Computer resources required by the numerous simulations shown in this paper were provided by the High Performance Computing Group at Los Alamos. Thanks are given to two anonymous reviewers whose comments led to significant improvements in the manuscript as well as to Jeremy Sauer, Vince Mousseau, Dana Knoll, and Piotr Smolarkiewicz for providing highly constructive comments over the past several years that helped in refining the ideas presented in this paper.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX

#### Microphysical Model

This appendix describes the chosen microphysical model that is composed of mass and number sources or sinks as well as a term that represents the falling of *ρ _{c}* and

*N*(see Table A1 for an overview of the components that compose the microphysical model). Note, however, that the intent of this paper is not to “sell” the current microphysical model. Rather, the purpose is to demonstrate that if a microphysical model is numerically differentiable and produces values that are physically realistic, then possible benefits of this effort are small temporal error growth during the course of a given simulation and the ability to reasonably agree with observed cloud fields.

_{c}The mass source *f*_{mass} in Eq. (5) is represented as follows:

In Eq. (A1), *f*_{qclimit} is a source term to keep *ρ _{c}* from falling significantly below the Newton tolerance and is expressed as follows:

where *ϕ*_{qclimit} = 1 s^{−1} is an inverse time scale and *q*_{cmax} ≈ 10^{−3} kg kg^{−1} is the expected maximum cloud water content within the computational domain with *f*_{qclimit} being nearly zero in value for *q _{c}* >

*q*. Likewise,

_{ce}*f*

_{Nclimit}is a source term to keep cloud number density

*N*or the number of activated cloud droplets per unit volume of air from falling below the background environmental value

_{c}*N*and is expressed in the same form as Eq. (A2) with the mass of the small cloud droplet being

_{ce}where *ρ _{w}* is the density of water with the radius of a new cloud droplet assumed to be 10

^{−6}m.

In Eq. (A1), Activate, represents a crude bulk parameterization (Twomey 1959) to convert aerosol number density, *N _{a}*, into cloud number density, and is expressed as follows:

where *ϕ*_{activate} = 1 s^{−1} is an inverse time scale for activation and *S*_{max} is the maximum supersaturation that defines the rate at which the target number of cloud droplets is reached. In Eq. (A4),

where *S*_{pos} is the positive part of the supersaturation field

and *S _{e}* = 0.1

*ε*

_{nl}is a value for which the supersaturation field can be reasonably resolved by the JFNK method.

The saturated mixing ratio *q _{υs}* found in Eq. (A6) is defined as follows:

with the saturated vapor pressure over water *e _{w}* expressed as

The term Deactivate in Eq. (A1) represents the rapid evaporation of *N _{c}* back into

*N*and is expressed, item 9 in Table 1 of Chen and Liu (2004), by the following relationship:

_{a}where

is the negative part of the supersaturation field, and *S*_{negsign} is used to smoothly turn off the parameterization for supersaturations greater than zero [i.e., neglect *S* after 0.5 in Eq. (A10)].

In Eq. (A1), Cond represents the diffusional growth of cloud droplets, item 1 in Table 1 of Chen and Liu (2004) obtained by statistical analysis of the condensational term from a bin model run in a parcel framework, and is expressed as follows:

where

is the inverse of the condensational time scale with the thermal conductivity of air, and *D _{υ}* = 0.015

*T*− 1.939 the diffusivity of water vapor.

For small amounts of *q _{c}* and

*S*< 0 the condensation model can remove more

*q*than is present within a given grid cell. To prevent this, a limiter has been implemented within the condensational model that has the following form:

_{c}where *ϕ*^{cond1} = *q*_{csmall}, and *ϕ*^{cond2} = 1 − *ϕ*^{cond1} with this limiter turning off evaporation when *q _{c}* <

*q*and

_{ce}In Eq. (6) *N _{c}* sources or sinks are represented as follows:

where the last term stands for the loss of *N _{c}* due to cloud droplet self-collection, item 3 in Table 1 of Chen and Liu (2004), and is expressed as follows:

with

and

is the mean cloud droplet radius with the size being limited as follows:

The parameterization for fall* _{ψ}* was computed by examining statistics from a Lagrangian parcel model where

represents the change of either *ρ _{c}* or

*N*due to the falling of cloud droplets with

_{c}Because the application of extra diffusion near cloud edges may lead to enhanced evaporation, the option exists in a simulation for evaporative processes through either Cond or Deactivate to be limited as follows:

where

with

and Deactivate being limited in the same manner as Cond. In Eq. (A23), *τ*_{evlimit} is a smooth representation for the ratio of the time scales *ϕ*^{ratio} produced by either the subgrid turbulent model, the cloud edge, or TKE with *ϕ*^{evtune} a tuning coefficient that incorporates the uncertainty associated with either time scale into the limiter.

## Footnotes

*Corresponding author address:* J. M. Reisner, Los Alamos National Laboratory, MS D401, Los Alamos, NM 87545. Email: reisner@lanl.gov