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.

Table 1.

Definition of key abbreviations used within the manuscript.

Definition of key abbreviations used within the manuscript.
Definition of key abbreviations used within the manuscript.

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

 
formula

where all indices range from 1 to 3; u1′ and u2′ are the Cartesian gas velocities in the horizontal, x1′ and x2′, directions; u3′ is the Cartesian gas velocity in the vertical, x3′, direction; Ω1′ = 0 with Ω2′ = 2Ω sinϕ and Ω3′ = 2Ωcosϕ being the x3′ and x2′ components of the earth’s rotation axis at the latitude ϕ; ρ is the gas density of the air, ρ = ρd + ρυ with ρd the dry air density and ρυ the water vapor density; p′ = ppe is the pressure perturbation with p the pressure of the gas computed using the ideal gas relationship [i.e., Eq. (11) in Reisner et al. 2005)] and pe = pe(z) the environmental pressure; ρ′ = ρρe is the density perturbation with ρe = ρe(z) the environmental density and the buoyancy force only being applied in the vertical direction, δi′3; ρc′ = qcρρce is the perturbation cloud water density with qc 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 ρce = qceρ, where qce = 0.1εnl qcb is the environmental cloud water mixing ratio with qcb = 1 kg kg−1 a constant base-state qc field; ueiare the environmental winds in all three directions; τij = (∂ui/∂xj) + (∂uj/∂xi) − (2/3)δijus/∂xs 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 ui are defined in Table 3.

Table 2.

Values of constants used in the cloud model.

Values of constants used in the cloud model.
Values of constants used in the cloud model.
Table 3.

Definition of key cloud model variables used within the manuscript.

Definition of key cloud model variables used within the manuscript.
Definition of key cloud model variables used within the manuscript.

The energy equation is expressed as

 
formula

where θ is the potential temperature θ = T(po/p)Rd/Cp, with T the temperature of the gas, fenergy = (θρL/TCp) fmass represents the release of energy during phase conversion of a given amount of mass fmass, associated with various microphysical source terms (described in the appendix); frad 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 x1′ direction is Fθ1′ = ρκ(∂θ/∂x1′).

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

 
formula

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.09LsTKE with Ls being the turbulence length scale.

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:

 
formula

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:

 
formula

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 ρc is similar to Eq. (5), aside from the term fmass being positive and the addition of a term to represent the fall of ρc.

The conservation equation for cloud number density Nc can be written as follows:

 
formula

where the fall term fallNcand the source term fnumber 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:

 
formula

where is a closure coefficient; the absolute value of a velocity component |ui| is approximated during a simulation by the following smooth function, |ui| ≈ ui tanh(100ui); Lextra is a length scale of the extra diffusion that is set to Δxi during a simulation; ε = 10−12 kg kg−1 is a constant to prevent division by zero; and Δiqc represents the change of the qc field in the xi direction. To ensure that the cloud-edge diffusion coefficient does not introduce fast time scales, the coefficient has been limited as follows:

 
formula

where

 
formula

with Δtref = 10Δtsound and Δtsound 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:

 
formula

and the second step is

 
formula

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 xi = iΔxi on the grid. Notice, as ϕCE−ADV increases in magnitude, tanh[ϕCE−ADV(ρc1/ρc*)] 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. Discrete model

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

 
formula

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:

 
formula

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:

 
formula

where εnl = 1 × 10−5 and the norm in Eq. (14) is defined by with ntot as the total number of relevant grid variables and as the iteration index.

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

 
formula

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υ, qc, Nc, 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 ϕCE−ADV = 1 for qυ, Nc, and TKE and ϕCE−ADV = 100 for qc; 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.

Table 4.

Key items in DYCOMS-II model setup of section 3.

Key items in DYCOMS-II model setup of section 3.
Key items in DYCOMS-II model setup of section 3.

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 Nc field that was initialized from data shown in Fig. 2 (Nc ≈ 100 cm−3) of Twohy et al. (2005). To initiate vertical motions, a weak moist bubble was introduced into the model domain of the following form:

 
formula

where ftime = (t − 330 s)/30 s; and fspace = dis(xi)/180 m with dis(xi) the distance in meters from the center of the bubble (x1′ = 1600 m, x2′ = 1600 m, and x3′ = 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.

Fig. 2.

Mean vertical profiles averaged over the last hour of (a) qc, (b) Nc, (c) variance of w, and (d) third moment of w from DYCOMS-II simulations SC1 (dashed redline), SC2 (red line), SC3 (green line), and TC1 (blue line). In (a) the dashed black line indicates both the initial and observed qc profile from flight 1 of DYCOMS-II.

Fig. 2.

Mean vertical profiles averaged over the last hour of (a) qc, (b) Nc, (c) variance of w, and (d) third moment of w from DYCOMS-II simulations SC1 (dashed redline), SC2 (red line), SC3 (green line), and TC1 (blue line). In (a) the dashed black line indicates both the initial and observed qc profile from flight 1 of DYCOMS-II.

A variable vertical resolution using 115 grid points produced by the following relationship (x3′ is in meters):

 
formula

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 qc, 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 qc and Nc (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.

Fig. 1.

Isosurfaces of qc = 10−4 kg kg−1 at the end of DYCOMS-II simulations (a) SC1, (b) SC2, (c) SC3, and (d) TC1.

Fig. 1.

Isosurfaces of qc = 10−4 kg kg−1 at the end of DYCOMS-II simulations (a) SC1, (b) SC2, (c) SC3, and (d) TC1.

Given that simulations SC3 and TC1 were of similar design, it is not surprising that their mean qc and Nc 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.

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 qc 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 qc at cloud top, but because of the falling term [e.g., Eq. (A20)], this greater amount of qc falls downward, thus producing a thicker cloud.

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 qc 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.

Fig. 3.

Mean vertical profiles averaged over the last hour of (a) qc, (b) Nc, (c) variance of w, and (d) third moment of w from DYCOMS-II simulations SC2 for ϕevtune = 10 (red line) and SC2 for ϕevtune = 1 (purple line). In (a) the dashed black line indicates both the initial and observed qc profile from flight 1 of DYCOMS-II.

Fig. 3.

Mean vertical profiles averaged over the last hour of (a) qc, (b) Nc, (c) variance of w, and (d) third moment of w from DYCOMS-II simulations SC2 for ϕevtune = 10 (red line) and SC2 for ϕevtune = 1 (purple line). In (a) the dashed black line indicates both the initial and observed qc profile from flight 1 of DYCOMS-II.

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:

 
formula

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.

Table 5.

Key items for the temporal error setup of section 4.

Key items for the temporal error setup of section 4.
Key items for the temporal error setup of section 4.

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Δtsound ≈ 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:

 
formula

was used to initialize a moist bubble where ftime = (t − 330 s)/30 s; and fspace = dis(x1′, x3′)/180 m with dis(x1′, x3′) the distance in meters from the center of the bubble (x1′ = 900 m, x3′ = 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 qc field) used to quantify differences between a given simulation and a reference simulation is defined as follows:

 
formula

with qcref coming from a reference simulation and qcl 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 qc 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 qc 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.

Fig. 4.

Two-dimensional plots of qc and wind vectors (m s−1) produced by the reference simulation, SC2–1/160r, revealing a moist bubble intercepting a stratus field at the following 60-s intervals: (a) 420, (b) 480, (c) 540, and (d) 600 s with shading for qc of 0.0001 (purple), 0.0002 (light blue), 0.0004 (dark blue), 0.0006 (orange), and 0.001 (green) kg kg−1.

Fig. 4.

Two-dimensional plots of qc and wind vectors (m s−1) produced by the reference simulation, SC2–1/160r, revealing a moist bubble intercepting a stratus field at the following 60-s intervals: (a) 420, (b) 480, (c) 540, and (d) 600 s with shading for qc of 0.0001 (purple), 0.0002 (light blue), 0.0004 (dark blue), 0.0006 (orange), and 0.001 (green) kg kg−1.

Fig. 5.

Two-dimensional plots of qc and wind vectors (m s−1) from (a) SC2–1/8 and (b) SC2–1/160 with the shading the same as in Fig. 4. (bottom) Two-dimensional plots of |qcqcref| and wind vectors (|uuref| and |wwref|) are shown for (c) SC2–1/8 and (d) SC2–1/160 with contour intervals of |qcqcref| of 1 × 10−6 (purple), 2 × 10−6 (light blue), 3 × 10−6 (dark blue), 6 × 10−6 (orange), and 10 × 10−6 (green) kg kg−1.

Fig. 5.

Two-dimensional plots of qc and wind vectors (m s−1) from (a) SC2–1/8 and (b) SC2–1/160 with the shading the same as in Fig. 4. (bottom) Two-dimensional plots of |qcqcref| and wind vectors (|uuref| and |wwref|) are shown for (c) SC2–1/8 and (d) SC2–1/160 with contour intervals of |qcqcref| of 1 × 10−6 (purple), 2 × 10−6 (light blue), 3 × 10−6 (dark blue), 6 × 10−6 (orange), and 10 × 10−6 (green) kg kg−1.

Fig. 6.

As in Fig. 5, but for TC1.

Fig. 6.

As in Fig. 5, but for TC1.

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).

Fig. 7.

Line plots of Error(qc) vs time from simulations (a) TC1–1/8 (black line), TC1–1/16 (red line), TC1–1/32 (green line), TC1–1/64 (blue line), and TC1–1/160 (orange line) discussed in section 4. (b) As in (a), but for SC2 with SC2–1/8 not being shown.

Fig. 7.

Line plots of Error(qc) vs time from simulations (a) TC1–1/8 (black line), TC1–1/16 (red line), TC1–1/32 (green line), TC1–1/64 (blue line), and TC1–1/160 (orange line) discussed in section 4. (b) As in (a), but for SC2 with SC2–1/8 not being 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.

Table 6.

Errors (×10−6) for qc computed at the end of the four simulation sets discussed in section 4.

Errors (×10−6) for qc computed at the end of the four simulation sets discussed in section 4.
Errors (×10−6) for qc computed at the end of the four simulation sets discussed in section 4.

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.

Table A1. Definition of key variables used within the microphysical model presented in the appendix.

Table A1. Definition of key variables used within the microphysical model presented in the appendix.
Table A1. Definition of key variables used within the microphysical model presented in the appendix.

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

REFERENCES
Andrejczuk
,
M.
,
J. M.
Reisner
,
B.
Henson
,
M. K.
Dubey
, and
C. A.
Jeffery
,
2008
:
The potential impacts of pollution on a nondrizzling stratus deck: Does aerosol number matter more than type?
J. Geophys. Res.
,
113
,
D19204
.
doi:10.1029/ 2007JD009445
.
Chen
,
J-P.
, and
S.
Liu
,
2004
:
Physically based two-moment bulkwater parameterization for warm-cloud microphysics.
Quart. J. Roy. Meteor. Soc.
,
130
,
51
78
.
Grabowski
,
W.
, and
P.
Smolarkiewicz
,
1990
:
Monotone finite-difference approximations to the advection-condensation problem.
Mon. Wea. Rev.
,
118
,
2082
2097
.
Grabowski
,
W.
, and
H.
Morrison
,
2008
:
Toward the mitagation of spurious cloud-edge supersaturation in cloud models.
Mon. Wea. Rev.
,
136
,
1224
1234
.
Jeffery
,
C.
, and
J.
Reisner
,
2006
:
A study of cloud mixing and evolution using PDF methods. Part I: Cloud front propagation and evaporation.
J. Atmos. Sci.
,
63
,
2848
2864
.
Klaassen
,
G.
, and
T.
Clark
,
1985
:
Dynamics of the cloud-environment interface and entrainment in small cumuli: Two dimensional simulations in the absence of ambient shear.
J. Atmos. Sci.
,
42
,
2621
2642
.
Leonard
,
B.
, and
J.
Drummond
,
1995
:
Why you should not use ‘hybrid’ ‘power-law’ or related exponential schemes for convective modeling—There are better alternatives.
Int. J. Numer. Methods Fluids
,
20
,
421
442
.
Margolin
,
L.
,
J.
Reisner
, and
P.
Smolarkiewicz
,
1997
:
Application of the volume-of-fluid method to the advection–condensation problem.
Mon. Wea. Rev.
,
125
,
2265
2273
.
Morrison
,
H.
, and
W.
Grabowski
,
2008
:
Modeling supersaturation and subgrid-scale mixing with two-moment bulk warm microphysics.
J. Atmos. Sci.
,
65
,
792
812
.
Mousseau
,
V.
,
D.
Knoll
, and
J.
Reisner
,
2002
:
An implicit nonlinearly-consistent method for the two-dimensional shallow-water equations with Coriolis force.
Mon. Wea. Rev.
,
130
,
2611
2625
.
Pielke
,
R.
,
1984
:
Mesoscale Meteorological Modeling.
Academic Press, 612 pp
.
Reisner
,
J.
,
S.
Wynne
,
L.
Margolin
, and
R.
Linn
,
2000
:
Coupled atmospheric–fire modeling employing the method of averages.
Mon. Wea. Rev.
,
128
,
3683
3691
.
Reisner
,
J.
,
A.
Wyszogrodzki
,
V.
Mousseau
, and
D.
Knoll
,
2003
:
An efficient physics-based pre-conditioner for the, fully implicit solution of small-scale thermally driven atmospheric flows.
J. Comput. Phys.
,
189
,
30
44
.
Reisner
,
J.
,
V.
Mousseau
,
A.
Wyszogrodzki
, and
D.
Knoll
,
2005
:
A fully implicit hurricane model with physics-based preconditioning.
Mon. Wea. Rev.
,
133
,
1003
1022
.
Skamarock
,
W.
,
2006
:
Positive-definite and monotonic limiters for unrestricted-time-step transport schemes.
Mon. Wea. Rev.
,
134
,
2241
2250
.
Smolarkiewicz
,
P.
, and
L.
Margolin
,
1998
:
MPDATA: A finite difference solver for geophysical flow.
J. Comput. Phys.
,
140
,
4459
4480
.
Stevens
,
B.
, and
Coauthors
,
2005
:
Evaluation of large-eddy simulations via observations of nocturnal marine stratocumulus.
Mon. Wea. Rev.
,
133
,
1443
1462
.
Stevens
,
B.
,
R.
Walko
,
G.
Feingold
, and
W.
Cotton
,
1996
:
The spurious production of cloud edge supersaturations by Eulerian models.
Mon. Wea. Rev.
,
124
,
1034
1041
.
Twohy
,
C.
,
M. D.
Petters
,
J. R.
Snider
,
B.
Stevens
,
W.
Tahnk
,
M.
Wetzel
,
L.
Russell
, and
F.
Burnet
,
2005
:
Evaluation of the aerosol indirect effect in marine stratocumulus clouds: Droplet number, size, liquid water path, and radiative impact.
J. Geophys. Res.
,
110
,
D08203
.
doi:10.1029/2004JD005116
.
Twomey
,
S.
,
1959
:
The nuclei of natural cloud formation. Part II: The supersaturation in natural clouds and the variation of cloud droplet concentrations.
Geofis. Pura Appl.
,
43
,
243
249
.
Zalesak
,
S.
,
1979
:
Fully multidimensional flux-corrected transport algorithm for fluids.
J. Comput. Phys.
,
31
,
335
362
.

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 Nc (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.

The mass source fmass in Eq. (5) is represented as follows:

 
formula

In Eq. (A1), fqclimit is a source term to keep ρc from falling significantly below the Newton tolerance and is expressed as follows:

 
formula

where ϕqclimit = 1 s−1 is an inverse time scale and qcmax ≈ 10−3 kg kg−1 is the expected maximum cloud water content within the computational domain with fqclimit being nearly zero in value for qc > qce. Likewise, fNclimit is a source term to keep cloud number density Nc or the number of activated cloud droplets per unit volume of air from falling below the background environmental value Nce and is expressed in the same form as Eq. (A2) with the mass of the small cloud droplet being

 
formula

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, Na, into cloud number density, and is expressed as follows:

 
formula

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

 
formula

where Spos is the positive part of the supersaturation field

 
formula

and Se = 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:

 
formula

with the saturated vapor pressure over water ew expressed as

 
formula

The term Deactivate in Eq. (A1) represents the rapid evaporation of Nc back into Na and is expressed, item 9 in Table 1 of Chen and Liu (2004), by the following relationship:

 
formula

where

 
formula

is the negative part of the supersaturation field, and Snegsign 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:

 
formula

where

 
formula

is the inverse of the condensational time scale with the thermal conductivity of air, and Dυ = 0.015T − 1.939 the diffusivity of water vapor.

For small amounts of qc and S < 0 the condensation model can remove more qc 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:

 
formula

where ϕcond1 = qcsmall, and ϕcond2 = 1 − ϕcond1 with this limiter turning off evaporation when qc < qce and

 
formula

In Eq. (6) Nc sources or sinks are represented as follows:

 
formula

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

 
formula

with

 
formula

and

 
formula

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

 
formula

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

 
formula

represents the change of either ρc or Nc due to the falling of cloud droplets with

 
formula

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:

 
formula

where

 
formula

with

 
formula

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