## Abstract

In this study, we propose a physical-process-based stochastic parameterization scheme using cellular automata for NOAA’s Next Generation Global Prediction System. The cellular automata, used to simulate stochastic processes such as the production and destruction of subgrid convective elements, are conditioned on unresolved vertical motion that follows a prescribed stochastically generated skewed distribution (SGS). The SGS is described by a stochastic differential equation and linked to observations by taking into account the first four moments from an observed dataset. In the proposed parameterization framework, we emphasize the need for a dynamical memory term to be included in physical-process-based stochastic parameterizations, and we illustrate the requirement for the dynamical memory using the Mori–Zwanzig formalism. Although this paper focuses on the methodology, early results indicate that if we apply our stochastic framework to deep cumulus convection, it is found that the frequency distribution of precipitation is improved in a single-member stochastic forecast, and some improved spread–skill relationship in ensemble runs can be found in state variables in the tropics, as well as in the subtropics.

## 1. Introduction

In numerical weather prediction (NWP), a set of nonlinear partial differential equations of fluid and thermodynamics are solved in order to predict the state of the atmosphere in the future, given an analysis of the present state of the atmosphere. Numerical methods are used to solve the equations on a discrete grid, and physical processes that are not resolved on the given grid or in time need to be represented by a parameterization. These approximations imply that even if a perfect initial condition could be provided to solve the equations, errors in weather and climate forecasts are inevitable. Over the past two decades, various stochastic parameterization schemes have been developed and applied at operational NWP centers around the world to address model uncertainty associated with model discretization and parameterization [e.g., Berner et al. (2017) and references therein].

The persistent interaction between resolved and unresolved dynamics and physics is often larger than can be represented by treating the stochastic term and the resolved dynamics as independent terms in a model. There are several ways of treating the effect of this lingering interaction, depending on the system at hand. One way is to stochastically approximate the dynamical system that allows memory of the subgrid dynamics on the resolved variables, a procedure referred to as “homogenization” (see, e.g., Gardiner 1985). In fact, as demonstrated by Gottwald et al. (2017), the homogenization technique can be expressed for the Mori–Zwanzig formalism, which has recently been regarded as useful for illustrative purposes because the dynamical memory and subgrid stochastic perturbations occur as separate terms (e.g., Evans and Morriss 2008; Franzke et al. 2015). A stochastic perturbation to the model-state tendency due to unresolved dynamics and physics at time *t* can be expressed under the Mori–Zwanzig formalism as the sum of two terms:

Two examples of this expression in a reduced form are the widely used stochastically perturbed physics tendencies (SPPT; Buizza et al. 1999; Palmer et al. 2009) and stochastic kinetic energy backscatter (SKEBS; Shutts 2005) schemes, developed at the European Centre for Medium-Range Weather Forecasts (ECMWF). In both these schemes, dynamical memory is approximated by a deterministic dissipation term, and subgrid stochastic perturbations are approximated by a Wiener process. In this paper, we develop a model framework that addresses the subgrid stochastic perturbation term through process-level perturbations on the model subgrid. However, the dynamical memory—or the large-scale spatial and temporal correlations of noise—is a crucial term in order to address the fact that local perturbations in physical processes have a memory in time and are influenced by scales far beyond the truncation of the model (e.g., Palmer 2001; Shutts 2013; Arnold et al. 2013; Ollinaho et al. 2017). Therefore, we cannot neglect this term in our development of a process-based stochastic parameterization, and we propose a framework in which we first perturb physical processes that contribute to the physics tendency and then perturb the physics tendency itself, using large-scale temporal and spatial correlations to address the fact that this state is dependent on the model state at earlier times.

We first present a discussion on dynamical memory, using the Mori–Zwanzig formalism to illustrate the essential terms in the stochastic parameterizations of uncertainties in subgrid dynamics. The aim is to make a connection between the theory that is well established in the literature of statistical physics and the stochastic physics parameterization development for NWP. In particular, we aim to argue why relatively large correlation scales are needed to address uncertainties associated with physical processes on the model’s subgrid.

### a. Dynamical memory

Stochastic models are often used to describe multiscale systems. In such models, resolved variables represent slow-system degrees of freedom, while noise represents the unresolved degrees of freedom that are assumed to vary rapidly. Often, the fast and slow variables in these models are represented as uncorrelated. However, such a formalism might neglect the interactions between fast and slow variables that can depend on the history, or dynamical memory, of this interaction. We illustrate the existence of dynamical memory below.

The Mori–Zwanzig formalism is well established in the literature of statistical physics and is applicable when external forcing does not depend on time (see, e.g., Evans and Morriss 2008). This formalism has recently been suggested as a theoretical foundation of stochastic parameterization development in climate modeling (see, e.g., Franzke et al. 2015). The Mori–Zwanzig projection operators project onto the linear subspace of slow-phase space functions, and the definition of such a projection allows separation of the Liouville equation by projecting onto the slow and the fast subspaces, which can be used to derive a reduced model in the form of a generalized Langevin equation. This is helpful in demonstrating how stochastic noise interacts with the resolved dynamics through a memory term. We will summarize the steps to get to the general Langevin equation below, while readers are referred to Evans and Morriss (2008), Givon et al. (2004), Gottwald et al. (2017), and Chorin and Hald (2013) for a more comprehensive derivation.

A model representing the nonlinear atmospheric system might be written as

where *t* is time, the dot denotes a time derivative, and **M** is the functional representation of the model. The model state is partitioned into the resolved () and unresolved () parts (i.e., ). Following Gottwald et al. (2017), one can rewrite Eq. (2) as the so-called Liouville equation, which is a linear partial differential equation. By first defining the Liouville operator (using tensor notation)

the Liouville equation associated with the model expressed by Eq. (2) can be used to describe an observable **z**(**x**, *t*) in the following form:

It can be shown that the solution of Eq. (3) is , where is the solution of Eq. (2) with initial condition . That is, if one can solve Eq. (3) for every , one can set = and obtain the *i*th component of for every *i*.

Using the Mori–Zwanzig projection operators with so-called semigroup notation (Chorin and Hald 2013), the solution of Eq. (3) can then be expressed as

where is a projector to map to the slow variables , and is a projector to map to the fast variables . The orthogonality between and is assumed for illustrative purposes, which is consistent with the practice of representing a model state as a point in a phase space to describe the model dynamics.

Following Givon et al. (2004) and Gottwald et al. (2017), one can then obtain

Equation (4) is exactly the same as Eq. (2), but in the form of the generalized Langevin equation.^{1} The first term in Eq. (4) is the dynamics of the resolved variables; the second term represents the “memory” term because it is an integration of quantities that are dependent on the model state at earlier times; and the third term is often called the “noise” term, representing the unresolved dynamics that is orthogonal to the resolved dynamics. Note that Eq. (4) applies to both deterministic and stochastic parameterizations. When applied to stochastic parameterizations, the sum of the last two terms in Eq. (4) represents the stochastic perturbation to the model-state tendency due to unresolved dynamics and physics that is introduced in Eq. (1).

The notion that there is an uncertainty at the subgrid scale, which consists of past uncertainty projected onto the current state, has also been discussed in the literature in terms of up- and downscale cascade of error growth. In the Lorenz (1969) idealized model framework for a system with a background kinetic energy spectrum that follows a *k*^{−5/3} power law, where *k* is the horizontal wavenumber, there is an insensitivity to the scale of the initial errors (Rotunno and Snyder 2008). As summarized by Weyn and Durran (2017, p. 2191), “This arises because large-scale errors propagate downscale very rapidly, which in turn saturate the error in the smallest scales, and then propagate back upscale as if they had simply been originated at the small scales.” In a global NWP model, however, such a *k*^{−5/3} power law is not realized with grid spacings of ~25 to ~50 km, as the mesoscale motions are not resolved. Thus, model error generated by subgrid-scale perturbations alone will not effectively propagate to the larger scales.

Furthermore, Durran and Gingrich (2014), Durran and Weyn (2016), and Weyn and Durran (2017) have shown that even when a *k*^{−5/3} power law exists in the KE spectrum, for certain weather regimes (i.e., U.S. snowstorm, deep moist convection, or mesoscale convective systems), it is difficult to identify such an upscale cascade of error growth. They state that the initial-condition perturbations that initially were strongest at the longest wavelengths “tended to grow simultaneously at all wavelengths, rather than through an upscale cascade, and that the error growth is rather ‘up amplitude’ than through a cascade from the smallest scales” (Weyn and Durran 2017, p. 2191); this is in line with the results we find in our present study.

In our framework, we do not simulate Eq. (4) exactly. However, we try as much as possible to implement the spirit of including each term. The first term is the resolved model. The last term, which represents the subgrid-scale effect on the large scale, is modeled as a combination of a cellular automaton and stochastic differential equations, the details of which will be provided below. The middle term represents the integrated effects of these subgrid-scale processes on the resolved dynamics. To some extent, this “memory” is included with the details of the cellular automaton setup, which allows intergrid information to propagate. However, the cross-boundary effects of the cellular automata are not sufficient to account for the integrated effect, so we supply an additional memory of subgrid effects by means of a first-order autoregressive process (AR(1)) perturbation to the tendencies. For the time being, we use the SPPT method to incorporate AR(1) in perturbations upon the tendencies from the physical parameterizations.

In SPPT, the model uncertainty due to physical parameterizations is represented by perturbing the net tendencies from the physical parameterizations of radiation, cloud physics, vertical diffusion, gravity wave drag, and convection by multiplicative noise as

Here, is the wind, temperature, or humidity of the numerical weather model; *μ* ∈ 0, 1 tapers the perturbations to zero near the surface and in the stratosphere; and is the first-order autoregressive AR(1) process for evolving spectral coefficients:

where controls the correlation over time step , and spatial correlations (Gaussian around the globe) for each wavenumber define for white random numbers *η*; *τ* is the temporal decorrelation scale.

Given the importance of the memory term as presented in this section, in our present development described in this paper, we adapt the perturbation on the tendencies following Eq. (5) from the SPPT scheme in order to include dynamical memory in our perturbations. The subgrid-scale perturbations are generated by cellular automata, which also to an extent exhibit decorrelation time scales through their self-organizational properties. In this development, we address mixing from deep cumulus convection only; however, we envision expanding the framework to also address uncertainties associated with physical processes such as turbulence, shallow convection, in-cloud microphysics, and radiation. This means that in our development, is the tendency due to convection only; thus, following the notation of Christensen et al. (2017), we call this independent (convection only) tendency perturbation “iSPPT (conv).”

### b. Subgrid stochastic perturbation

The most common method to parameterize cumulus convection in today’s weather and climate models is based on the mass flux concept that was first formulated by Ooyama (1971), Yanai et al. (1973), and Arakawa and Schubert (1974). In the bulk mass flux concept, a single steady-state updraft is assumed to represent the effects of all active cloud elements in a grid box, averaged over the different stages of their life cycle. As the horizontal resolution of the NWP model increases, however, the finite number of convective updrafts within a grid box will generally be reduced. At ~10–15-km horizontal resolution, there may be about one active cloud element on average, although most realizations will very likely contain more or fewer cloud elements than this number. The ensemble-averaged effects of the subgrid cloud field may thus be very different from those of a particular realization; hence, the subgrid-scale variability should be sampled, rather than represented by an ensemble mean (Plant et al. 2015; Monahan and Culina 2011). To this end, a number of studies have explored whether stochastic approaches for deep convection can be used to mimic statistical fluctuation in cloud numbers and intensities (e.g., Lin and Neelin 2002; Bright and Mullen 2002; Grell and Dévényi 2002; Byun and Hong 2007; Plant and Craig 2008; Khouider et al. 2010; Tompkins and Berner 2008; Dorrestijn et al. 2013; Gottwald et al. 2016; Frenkel et al. 2013; Bengtsson et al. 2013; Sakradzija et al. 2016; Shutts 2015).

Craig and Cohen (2006a) used a statistical mechanics method to investigate the equilibrium fluctuations of convective cloud properties (such as mass flux) around an average taken over a limited region of convective clouds. It was hypothesized that the distribution of cloud number within a given area should follow that of a Poisson distribution, an idea that was supported through analysis of equilibrium cloud-resolving model simulations (Craig and Cohen 2006a; Cohen and Craig 2006b; Plant and Craig 2008; Plant et al. 2015). In recent years, stochastic birth/death models for cloud population have been explored within the NWP community. Such models are useful for simulating competing processes seen in nature, such as production versus dissipation, and the distribution of the population is Poisson-like. Plant (2012) proposed an idealized stochastic model framework of the cloud population in a grid box using a so-called master equation to evolve the probability distribution of convective cloud properties. In that study, transition probabilities representing births and deaths as well as environmental destabilization and stabilization were used to evolve the master equation. In a similar framework, Hagos et al. (2018) also considered a probabilistic representation for the nonequilibrium dynamics of cloud populations to predict the spectrum of cloud sizes and their evolution due to a given forcing.

Khouider et al. (2010), Dorrestijn et al. (2013, 2015, 2016), Gottwald et al. (2016), and Frenkel et al. (2013) also use birth/death processes to describe the population of clouds within a finite area. In these studies, Markov chains were used, and transition probabilities were given to model how a cloud can transition between different cloud types within a model grid box. In Dorrestijn et al. (2016), the population of these cloud types averaged over the area (model grid box) was then used to provide the convective area fraction for computing the convective mass flux at cloud base in the SPEEDY model. Closely related is the work by Bengtsson et al. (2011, 2013) and Bengtsson and Körnich (2016), where birth/death processes from cellular automata were used on the subgrid of the NWP model grid. The self-organizational nature of the cellular automaton allows for subgrid elements to form clusters not only on the subgrid, but also across the NWP model grid boundaries, allowing for communication of the cellular automaton across adjacent NWP grid boxes. In Bengtsson et al. (2013), the convective area fraction is given by the fraction of the cellular automaton cells that are “alive” within a grid box, which is further used in the prognostic equation of area fraction in the “3MT” deep convection scheme of the European regional ALARO model (Gerard et al. 2009; Termonia et al. 2018). Together with a prognostic equation for the updraft vertical velocity, the updraft area fraction modulated by the cellular automaton yields the mass flux.

In this study, we take advantage of the self-organizational properties of cellular automaton described by Bengtsson et al. (2013) to represent the subgrid stochastic term in Eq. (1). We divide the NWP model’s grid box into several subgrid elements and let the cellular automata act on this finer grid, such that the cellular automaton cells represent subgrid convective cloud elements. Physically, the existence of competing birth/death processes in the cellular automaton model is critical to make the distribution of the convective cloud numbers a skewed distribution.^{2} We use this distribution to provide, in a stochastic manner, the number of subgrid convective plumes in the deep convection parameterization of the model. An advantage of the cellular automaton model as implemented is that through its self-organization, subgrid cells can be triggered in an adjacent NWP model grid box, allowing for horizontal communication of the deep convection scheme, which is otherwise column based.

The transition rules of the cellular automaton are conditioned on a stochastically perturbed large-scale state, where the perturbations are given by a general class of non-Gaussian distributions called stochastically generated skewed (SGS; Sardeshmukh et al. 2015). The SGS probability density function can be associated with a stochastically perturbed damped linear Markov process and formulated as a stochastic differential equation. Thus, by using an SGS distribution for the perturbations of the large-scale model field, we do not have to sample from an existing PDF, but can rather use the stochastic differential equation directly for the perturbations. Furthermore, the distribution provided by the stochastic differential equation is made to be realistic by prescribing the first four moments of the distribution from high-resolution LES simulations or observational data.

We refer to this approach as CA-SGS, where the idea is to provide a physically based stochastic estimate of the subgrid variability of convection and introduce subgrid-scale organization. We explore the idea of using CA-SGS in one of the cumulus convection schemes developed for NOAA/NCEP’s next-generation Global Forecast System (GFS), as described in the next two sections.

## 2. Cellular automata and stochastically generated skewed distribution (CA-SGS)

A cellular automaton is a simple mathematical model useful for simulating the behavior of complex physical systems exhibiting competing processes, such as production versus dissipation of cumulus convection. In the cellular automata, space and time are discrete, and the state of a cell at the next time step depends on the state of its neighborhood at the current time step. Transition rules to evolve a cell from one state to another, based on the state of the neighborhood, are prescribed. This transition rules can be given either as a Boolean condition or as real numbers representing the probability for a cell to have a given state [lattice Boltzmann method (LBM)]. The latter has been shown to be successful for modeling high Reynolds number flows (Chopard and Droz 1998).

As elaborated on by Berner et al. (2017), applications of cellular automata for NWP have previously been explored within the atmospheric science community. In Berner et al. (2008), cellular automata were used as a stochastic pattern generator for SKEBS based on an initial idea presented by Palmer (2001) and, as described above, by Bengtsson et al. (2013) and Bengtsson and Körnich (2016) for use as a stochastic parameterization of cumulus convection. In all of these studies, the updating rules of the cellular automaton are referred to as “generations,” which is an extension to the famous cellular automaton “game of life” (proposed first by J. Conway in 1970); the intent is that the spatial and temporal scales and the continuous self-organization mimic the scales and the behavior observed from organized deep convection in the atmosphere.

In this study, we stay within the “generations” family of cellular automata in that we use a deterministic ruleset for the “birth” and “survival” of a cell on the subgrid of the NWP model grid following Conway’s game of life. However, an important difference from Bengtsson et al. (2013) and Bengtsson and Körnich (2016) is that the ruleset is conditioned on a stochastically perturbed large-scale state, described below. The reason why we choose to condition the cellular automata on a perturbed large-scale state rather than the gridbox mean of this state is to capture the subgrid variability of unresolved motion and to be able to inform our development of stochastic physics from observations. In addition, this transition condition for the cellular automaton adds a stochastic component in its evolution, compared with a strictly Boolean ruleset. The condition for our cellular automaton consists of subgrid perturbations of level-mean large-scale vertical velocity and/or specific humidity , with perturbations sampled from an SGS probability density function (PDF; Sardeshmukh et al. 2015). At the initial time, an active cellular automaton cell on the subgrid is then defined as the perturbed vertical velocity larger than a given threshold. Here, is the gridbox mean vertically averaged velocity between the surface and 350 hPa from our NWP model. This large-scale state was chosen following Dorrestijn et al. (2016) because it is known to have a very high correlation with deep convection activity (Peters et al. 2013; Davies et al. 2013)

In this study, we modify the SGS distribution to represent the interaction between the dynamical noise with the mean state and its effect on the anomalies. In the example of gridbox mean of the level mean vertical velocity, we add from the NWP model into Eq. (10) of Sardeshmukh et al. (2015), which now takes a slightly modified form:

Here, *E*, *g*, and *b* are the parameters of the distribution, and are random Gaussian variables with zero mean and unit variance that are uncorrelated both spatially and temporally, is the vertical velocity kept from the previous time step, and is a damping time constant viewed as a tuning parameter of the system. The difference from Sardeshmukh et al. (2015) is that the constant *g* has been replaced by to ensure that the subgrid vertical velocities are perturbed around the current gridbox mean. Following Sardeshmukh et al. (2015), the parameters *E*, *g*, and *b* can be determined from the first four moments of an observed sample series of At the first time step, we assume = 0 and simply set = as a crude approximation ( is the variance). Numerically, the equation is solved using a predictor-corrector method (Rümelin 1982).

In our study, we use the first four moments of vertical velocity from an observational dataset at Darwin, Australia, to compute the parameters *E*, *g*, and *b*. Future studies will include LES to compute the moments of vertical velocity to assess the robustness of these observations. A value from the distribution is sampled in each grid cell of the higher-resolution cellular automaton grid, using different random values of and in each cell.

In the case of level mean vertical velocity, is perturbed within each cellular automaton grid cell following

and a convective element is defined where . We do the same for specific humidity and explore the effects herein of conditioning the cellular automaton on a combination of these large-scale fields.

Figure 1 shows a schematic image of the cellular automaton grid (thin blue lines) overlaying the NWP model grid (thick gray lines). Each dark gray box represents a subgrid convective element (defined as ), where the state of this element at the next time step is determined using Conway’s game of life rules modified to include conditioning on each time step. Our cellular automaton is specifically described as follows:

If the state of the cell is 0, a new cell is born if or if the cell has exactly two or three neighbors meeting this criterion.

If the state of a cell is >1, it will survive if it has exactly four or five neighbors with the criterion . Otherwise, the state will be reduced by 1 from its assigned lifetime described below.

If a cell is born, it gets assigned a lifetime, which is an integer value equal to , where is a random number between 0 and 1, and nlives is the maximum prescribed lifetime a cell can have.

Normally, in deep convection parameterizations, horizontal communication only takes place via gridscale circulations. However, in nature, small-scale processes such as gravity waves and cold pool dynamics can act to trigger and organize convection in a way not possible by a one-dimensional plume model under a steady-state assumption (Huang 1990; Mapes and Neale 2011). The cellular automaton provides a way to communicate across neighboring NWP gridbox boundaries through its self-organizational behavior. This is possible since the cellular automaton grid acts on the subgrid of the numerical model (Bengtsson et al. 2013). Other examples of parameterization aimed to improve convective organization are those of Mapes and Neale (2011), where rain evaporation was used as a way of addressing convective memory, and Grandpeix et al. (2010), who proposed a density current parameterization to describe cold pools.

The distribution of convective elements on the subgrid depends on the chosen number of subgrid cells, the threshold criteria, and the transition rules of the cellular automaton (including which variable the cellular automaton is conditioned on and the choice of “lifetime” parameter). As a relevant example, we consider the case where our CA-SGS is acting on subgrid mesh with 10 × 10 points per NWP model grid box and a threshold criterion of 0 m s^{−1}. In Fig. 2a, the distribution of the subgrid perturbations from the SGS is shown on the finer cellular automata grid, and in Fig. 2b, the distribution of the subgrid number of convective elements provided by our model is shown. The numbers are given on the grid of the NWP model, and both distributions are heavily skewed.

Furthermore, the evolution of the cellular automaton pattern depends strongly on whether we condition the transition from one state to another on a model field such as vertical velocity, as can be seen in Fig. 3. In the left panel is the cellular automaton averaged back onto the NWP model grid after running freely, given an initial condition of convective elements defined by our SGS perturbed vertical velocity field. The initial condition is lost after only a few time steps, and the pattern becomes global and uniform in space. In the right panel is the cellular automaton conditioned on perturbed vertical velocity used to associate the pattern evolution with the activity of convection, as given by the ruleset described above. Clearly, the pattern is dominated by areas of positive large-scale vertical velocity; however, within the large-scale pattern, smaller-scale variability in the pattern is present.

## 3. Coupling of the CA-SGS to the Chikira–Sugiyama deep convection parameterization

NOAA GFDL’s finite-volume cubed-sphere (FV3) dynamical core has been selected for the new Next Generation Global Prediction System (NGGPS) atmospheric model, and this dynamical core using updated GFS model physics, hereafter referred to as FV3GFS, will replace NOAA’s current GFS. One of the cumulus parameterizations considered for the NGGPS is the Chikira–Sugiyama convection scheme (Chikira and Sugiyama 2010). In this scheme, the entrainment rate varies vertically depending on the surrounding environment. An entraining plume model is adopted based on Gregory (2001). Cloud types of different depths are spectrally represented, following the spirit of the Arakawa–Schubert scheme (Arakawa and Schubert 1974), according to different values of updraft velocity at cloud base. Furthermore, the evolution of the cloud-base mass flux for each plume type is predicted (rather than being diagnosed) using a prognostic convective kinetic energy closure (Randall and Pan 1993; Pan and Randall 1998).

As part of its implementation into the FV3GFS model, the scheme has been made scale adaptive, given the framework of Arakawa and Wu (2013), and has been tested together with the simplified higher-order closure (SHOC) planetary boundary layer (PBL) scheme described by Bogenschutz and Krueger (2013) and the two-moment cloud microphysics scheme described by Morrison and Gettelman (2008). This is the configuration of model physics used in this study.

The Chikira–Sugiyama scheme, like any other parameterization of subgrid processes, approximates the space–time average effect of many possible realizations of the subgrid moist convection. While the scheme is scale aware in that the space–time average effect of subgrid moist convection diminishes as the grid size decreases to subkilometer scales, the space–time average approximation cannot represent the average effect of natural fluctuations of subgrid convection. Therefore, a stochastic subgrid parameterization is needed to account for the stochastic effect of the fluctuations.

In the present implementation of the Chikira–Sugiyama scheme in the FV3GFS model, a spectrum of subgrid plumes is represented by different values of vertical velocity at cloud base. The authors of the scheme argue that cloud base properties influence the vertical profiles of in-cloud properties, as well as cloud-top height (a weaker upward motion of cloud air at cloud base leads to larger entrainment rate and smaller in-cloud moist static energy, thereby leading to a lower cloud top (Chikira and Sugiyama 2010, p. 67). The spectrum of cloud-base updraft velocities is assumed to first be given from the minimum to the maximum with a fixed interval, such as

where *N* is the number of plumes, and are set to fixed values to set the bounds of the range in the scheme, and *N*_{max} is the maximum number of plumes.

However, not all these possible plumes are generally considered to be active by the scheme—the actual number depends on the degree of convective instability in the environmental sounding.

Thus, if we want the CA-SGS model to provide the number of plumes, we cannot simply replace the plumes summed over in the Chikira–Sugiyama scheme, as this depends on environmental conditions. Therefore, we use the CA-SGS model to change the number of plumes by perturbing the environment (triggering criteria), as well as the maximum number of subcloud elements *N*_{max} used, which changes the interval of plumes used in Eq. (9).

In future developments, we plan to take advantage of the subgrid cloud size distribution that is effectively given by the CA-SGS, associating each subgrid cluster of convective elements with a particular cloud size or radius. However, since the current spectral representation of plume types in the current Chikira–Sugiyama convection scheme does not consider cloud radius to influence the entrainment, we simply (as a first step) let the number of active subgrid cloud elements from our CA-SGS model alter the number of maximum cloud elements *N*_{max} and then use it as a pattern to perturb the environment. The self-organizational aspects of the cellular automata are still of value, as we will sum over more cloud elements (and yield stronger updrafts) in areas where there is more organization.

The Chikira–Sugiyama scheme does not have an explicit triggering criterion; the convection is active in regions where the updraft vertical velocity, adapted from Gregory (2001), is positive:

where and *B* are the updraft velocity and the buoyancy of the cloud air parcel, respectively, and and are dimensionless constant parameters ranging from 0 to 1.

To perturb the environment, we adopt a method similar to that proposed by Tompkins and Berner (2008) to perturb parcel humidity in order to represent subgrid-scale thermodynamic variability not accounted for by using gridbox average quantities. Perturbing the boundary layer humidity (and/or temperature), and thereby the parcel buoyancy, will influence both the triggering and the strength of convection. Furthermore, as discussed by, for example, Tompkins and Berner (2008), Parsons et al. (2000), and Lopez and Moreau (2005), boundary layer humidity perturbations have greater effect than temperature perturbations. Thus, for simplicity, we only consider humidity perturbations. Other stochastic parameterizations addressing convective triggering include modifications to the Kain and Fritsch cumulus parameterization (Kain 2004), by stochastically perturbing the triggering condition, as in Bright and Mullen (2002) and Song et al. (2007). Another method was explored by Rochetin et al. (2014) using LES to analyze stochastic aspects of the convection initiation in the Emanuel (1991) parameterization.

In our scheme, instead of randomly sampling subgrid humidity distributions provided by the cloud scheme to perturb the parcel’s humidity (as in Tompkins and Berner 2008), we utilize the CA-SGS distribution to perturb the parcel humidity. For our perturbations to the humidity field, we condition the cellular automaton on vertical velocity and specific humidity in order to have a pattern different from that used to perturb the number of convective plumes. We use the pattern given by averaging the cellular automaton cells back onto the NWP model grid; we reduce the distribution to have a zero mean and then scale the amplitude such that the distribution ranges between −0.001 and 0.001. We finally perturb the input moisture to the deep convection scheme as

here is the gridbox mean specific humidity, and is the perturbation term. For each grid box, the same perturbation is used throughout the atmospheric column.

## 4. Sensitivity of the CA-SGS parameters

While there are many different matrices one can use to demonstrate the sensitivity to the parametric choices and the coupling strategies used in this study, we here look at the sensitivity of some of the CA-SGS parameters using a single case study in terms of ensemble spread and skill. A more complete evaluation of the stochastic scheme will be presented in section 6. Figure 4 shows the tropical domain average 200-hPa zonal wind (left) and 850-hPa temperature (right) root-mean-square error (RMSE) and ensemble spread of a 20-member ensemble using FV3GFS at C96 (~100 km) resolution for the forecast initialization 1 August 2014. The RMSE is computed comparing the model output to the ECMWF interim reanalysis (ERA-Interim; Dee et al. 2011). The solid lines represent the RMSE, and the dashed lines represent the ensemble spread (standard deviation). The different colors show the impact on ensemble spread and ensemble skill given different parametric choices or implementation/coupling strategies. All of the members have initial condition perturbations. The black curve does not have any stochastic physics switched on and is thus the reference experiment (labeled “control”). The red curve is the CA-SGS experiment that we have outlined so far in this paper (labeled “CA_control”), with ncells = 5 and nlives = 30 (model time step = 900 s), and all the other lines in Fig. 4 represent *single* alterations to this CA-SGS control experiment. [Note that in these sensitivity experiments, we do not apply AR(1) perturbations onto the deep convective tendencies; this sensitivity will be examined more carefully in section 6.]

First, we look at the sensitivity of spread and skill to the transition rules of the CA-SGS to see the impact if we do not condition the cellular automaton on the SGS-perturbed large-scale fields. In this case, the cellular automaton strictly follows the ruleset of Conway’s game of life (green curve labeled “CA_global”). The cellular automaton in this case looks like the pattern shown in the left panel of Fig. 3. The impact using this global pattern on ensemble spread is very small, much smaller compared to the experiment where we condition the cellular automaton on a large-scale field. One reason may be that conditioning on a large-scale state yields much larger correlation space and time scales in the perturbation pattern. It should be noted, however, that the scales of the cellular automaton without a large-scale condition can be larger than what is shown in Fig. 3. This can be achieved by, for instance, decreasing the number of subgrid cells or even running the cellular automaton on the native NWP model grid, as well as by increasing the lifetime parameter nlives. However, for a comparison with the same subgrid cells (ncells) and same nlives, it is clear that the effect on ensemble spread is greater if we condition the cellular automaton on a perturbed large-scale field. Furthermore, the cellular automaton conditioned on the SGS-perturbed vertical velocity shows an enhanced realism of convective regions. In particular, the right side of Fig. 3 shows clear signatures of orographic forcing and convergence zones. This realism is not obvious from global measures of skill.

After this, we look at the sensitivity to only perturbing the closure by addressing the convective plumes and remove the trigger perturbations to the input specific humidity following Eq. (11). This sensitivity is shown by the blue curve labeled “CA_closure_only.” As can be seen, the trigger perturbations are important for ensemble spread, since we otherwise only perturb the closure (strength) of the convection and do not allow any new convection to be triggered in neighboring grid cells.

Last, we investigate the impact of using a skewed distribution when perturbing our large-scale fields used to condition the cellular automaton. In this case, we replace the SGS distribution by a Gaussian distribution. However, it is important to note that even though the distribution with which we perturb the model large-scale field is now given by a Gaussian white noise, the cellular automaton “plume” number distribution will still be skewed because of its memory. Thus, it is not straightforward to assess the implications of this change by simply looking at ensemble/spread skill. Nevertheless, the experiment with replacing the SGS-distribution with a Gaussian distribution is given by the cyan-colored line labeled “CA_gaussian,” and as can be seen, there is very little impact on ensemble spread and skill using this change. We believe there is still a value in using the SGS distribution for our large-scale perturbations, as other studies have shown that the subgrid distribution of convective plumes more resembles a skewed distribution such as SGS or Poisson than a Gaussian (Craig and Cohen 2006a; Cohen and Craig 2006b; Plant and Craig 2008; Plant et al. 2015), and it gives an opportunity to inform our development of stochastic physics from observations.

## 5. Experimental design

We evaluate the impact of the proposed scheme in FV3GFS by generating both single-member stochastic and ensemble weather forecasts. As discussed in section 1, we argue that the large-scale memory term of Eq. (1) is important, and in this study, we use the first-order autoregressive process AR(1) from the SPPT scheme implemented in the FV3GFS model to represent this term by perturbing the tendencies given by the revised (stochastic) implementation of the Chikira–Sugiyama scheme outlined above. Table 1 summarizes the setup of the experiments.

The impact of the scheme is evaluated for 5-day weather forecasts for both single-member and ensemble integrations. Initial conditions were generated using the original spectral GFS model, as opposed to running a full data assimilation system. The ENS_CONTROL run does not have any model error representation using stochastic physics but uses initial condition perturbations.

The CA-SGS runs were performed using a vertical-velocity-conditioned cellular automaton with 5 × 5 subgrid cells and a threshold of 0 m s^{−1} for activation of convective elements. The cellular automaton is spun up by performing 1000 iterations at the first time step, where the maximum prescribed lifetime of one cell is (nlives × *dt*), where nlives = 30 and *dt =* 450 s*.*

We assess the impact of the scheme by comparing reforecast temperature, geopotential height, and winds to ERA-Interim (Dee et al. 2011). The predicted precipitation is compared against the Tropical Rainfall Measuring Mission (TRMM) 3B42 high-resolution data (3 hourly, 0.25° spatial resolution; Huffman and Bolvin 2011). As commonly seen in satellite precipitation estimates, TRMM3B42 has been determined to have large relative errors at small precipitation rates; however, time/area averaging significantly reduces the random error (Huffman et al. 2007; Huffman and Bolvin 2011).

## 6. Impact on model performance

Because of the space–time filter in NWP models, large to extreme amounts of precipitation are often underestimated compared with instantaneous point observations, particularly when the horizontal grid spacing is large, and thus vertical motions are not well represented. Since the stochastic parameterization presented here is developed to implicitly represent the variability of subgrid-scale fluctuations of deep convection, we first look at the frequency distribution of precipitation in a single-member stochastic run on weather forecast time scales to see if variability of precipitation is better represented. This same sort of assessment was recently conducted by Watson et al. (2017), who looked at the ECMWF stochastic physics developments SPPT and SKEBS, evaluating the tropical rainfall variability on weather forecast time scales in the Integrated Forecast System (IFS). Meanwhile, Wang and Zhang (2016) and Wang et al. (2017) looked at how using the Plant–Craig stochastic convection scheme (Plant and Craig 2008) affected the frequency distribution of tropical precipitation in the NCAR CAM5 model for various resolutions and time scales (both weather forecasts and climate simulations). In all these studies, it was found that using a stochastic parameterization improves the frequency distribution of tropical precipitation, in particular for large to extreme amounts.

Figure 5a compares the normalized frequency distribution of TRMM rainfall observations with the unperturbed control run (CONTROL), a single-member run using iSPPT on deep convection (iSPPT conv), and a single-member run using CA-SGS described above (iSPPT conv + CA-SGS). Both the iSPPT on convection and the CA-SGS stochastic physics increase the precipitation variability across the bin distributions. This increase is larger using CA-SGS, which may be explained by the fact that we now perturb processes within the physics with the cellular automaton pattern, as well as the physics tendencies that are output of the scheme with the AR(1) pattern. Furthermore, treating processes within the convection scheme in a stochastic manner allows for new convection to be triggered, which is not the case when using iSPPT alone: if the tendency is zero, the perturbation of the tendency will also be zero. In addition, perturbing the mass flux within the convection scheme (which is essentially what is done using CA-SGS to perturb the cloud number in each grid box) allows for a more direct perturbation to the convective precipitation, compared with iSPPT where only the temperature, humidity, and wind tendencies out of the convection scheme are perturbed. Similar reasoning was put forth by Ollinaho et al. (2017) and Leutbecher et al. (2017), where parameters within the physics were perturbed using the stochastic perturbed parameter (SPP) scheme in the ECMWF IFS model. The difference here is that we have attempted to address processes rather than parameters. Another interesting impact on the precipitation variability is that the problem of too much drizzle is also improved with CA-SGS by reducing the smallest threshold amounts compared with CONTROL and iSPPT.

In Fig. 5b, the standard error and mean bias are shown for the same experiment as described above. There is a small increase of the bias compared with CONTROL using the stochastic schemes, and for longer lead times, this bias is similar in both the iSPPT experiment and the iSPPT+CA-SGS experiment. The error is slightly increased using iSPPT and increased further using iSPPT+CA-SGS. Thus, it appears that the CA-SGS increases the error due to the increased variability of the predicted precipitation. The bias using iSPPT may come from the fact that only perturbing convection tendencies without perturbing the remaining physics tendencies may cause inconsistencies that introduce biases; this should be investigated further in the future. The improvement seen in Fig. 5a suggests that the standard error is not the best way of evaluating precipitation skill, as there can be a double penalty due to location error in a more variable forecast of large amounts. A careful assessment of precipitation skill against local observational networks needs to be carried out in the future using metrics such as fraction skill score, frequency bias, probability of detection, and probability of false alarm rate. The increase in bias is a concern, and it seems mainly (but not solely) related to the iSPPT perturbations, which need to be addressed in future work.

The stochastic nature of the cumulus parameterization will lead to a different response in the resolved scale in each ensemble member. Therefore, it is of interest to investigate the impact the scheme has on ensemble spread. A commonly used metric to understand if an ensemble system is over- or underdispersive is the spread–skill relationship. A perfect ensemble in this regard will have a one-to-one ratio between the standard error of the ensemble mean and the standard deviation of the ensemble spread (e.g., Whitaker and Loughe 1998; Bengtsson et al. 2008).

Figure 6 shows the 5-day forecast ensemble spread and error for the vertical profiles of temperature and zonal wind averaged over global longitude bands. The top panel shows the difference between ENS_iSPPT and the unperturbed ensemble control. The ensemble spread is increased, and the ensemble mean error is reduced in ENS_iSPPT, in particular over the tropics in the zonal wind field. When we turn on the stochastic deep convection parameterization through CA-SGS in addition to iSPPT, the ensemble spread is increased further, and the ensemble error is reduced over the tropics, with the exception of a degradation in error at about 100 hPa. Over the subtropics, the ensemble spread is increased compared to iSPPT with a mostly neutral impact on skill; however, there is a tendency of a degradation in error toward the poles. The CA-SGS without AR(1) perturbations onto the deep convective tendencies compared to the control is shown in the bottom panel. The CA_SGS + iSPPT experiment (middle panel) is not exactly the sum of the individual iSPPT (top panel) and CA-SGS (bottom panel) experiments since new convection can be triggered from the moisture perturbations by the CA-SGS.

The results from Fig. 6 can be viewed in terms of single-level domain averages in order to better visualize how the ensemble spread is related to the ensemble error. Figure 7 shows the 850-hPa zonal wind and temperature ensemble spread and skill (standard error) as a function of lead time for the ensemble experiments listed in Table 1. The main impact of the stochastic physics schemes is, not surprisingly, over the tropics. The dynamical memory is important to generate large ensemble spread, as indicated by the iSPPT experiment, and in addition to iSPPT the CA-SGS can add a further increase in ensemble spread and a small reduction in the ensemble mean error.

Next, we look at the model bias comparing our experiment outputs to ERA-Interim. The model temperature bias at forecast lead time 120 h (day 5) is shown in Fig. 8 as a vertical distribution of the zonal mean. The CA-SGS experiment tends to warm the atmosphere, which leads to a reduction of an existing cold bias in the control experiment in the 400–150-hPa layer (upper-left panel). However, it also leads to an amplification of an already-existing warm bias around 500 hPa. The iSPPT experiment has no impact on the model temperature bias. The uniform increase of temperatures using the CA-SGS scheme may come from the skewed distribution from the cellular automata used to perturb the input specific humidity, and this warrants further investigation.

Last, the impact on the spectral representation of kinetic energy ensemble spread and error is studied in Fig. 9 to investigate the uncertainty as a function of atmospheric scales. The left panel shows all the wavenumbers, whereas the right panel is zooming in on the tail of the spectra. Comparing different lead times (not shown), there is no evidence in an upscale or downscale cascade of model uncertainty, but rather, as discussed in section 1, there is an up-amplitude growth of the uncertainty at all wavenumbers simultaneously, which is in agreement with other studies (Durran and Gingrich 2014; Durran and Weyn 2016; Weyn and Durran 2017). The ENS_CA-SGS [without AR(1) perturbations] shows a small increase in the KE spectral ensemble spread at all wavenumbers, even though the scales on which the perturbations are inserted are close to the truncation scales of the model. The AR(1) perturbations, ENS_iSPPT (associated with large space/time scales correlations), are dominating the increase in ensemble spread at the largest scales, whereas the increase in ensemble spread at the smallest scales are the same for ENS_iSPPT and ENS_CA-SGS. The combination of AR(1) perturbations and CA-SGS (ENS_CA-SGS+iSPPT) increases the ensemble spread at all scales, but most efficiently at the smallest scales (in a relative sense, as the magnitude of the increase is greater at the larger scales). One reason for this could be that there is a downscale cascade of uncertainty given by iSPPT in combination with the local perturbations.

## 7. Conclusions and outlook

In this study, we proposed a physical-process-based stochastic parameterization using cellular automata in a model framework. In this framework of parameterization, we first perturb physical processes acting on the model subgrid scales. Then, the tendencies themselves are perturbed using a random spatial pattern with larger decorrelation scales to address the fact that the uncertainty is dependent on the model state at earlier times and at a distance from where the process is occurring.

The question is often raised as to why stochastic parameterizations need synoptic-scale decorrelation time and space scales to adequately address uncertainties associated with physical processes occurring at the model’s truncation scale. One common explanation is that some amount of added or implicit numerical diffusion is typically needed, such that local perturbations near the model truncation are rapidly damped. However, in practice, at operational centers, decorrelation scales well beyond the model’s “effective” resolution are chosen in order to achieve the best results.

One point that we have emphasized is that there is a need for a dynamical memory term in *any* sort of physical parameterization. Up to now, the importance of this memory term, which is often discussed in statistical physics, has been given little emphasis in the discussion of stochastic parameterizations implemented at operational centers. There is memory in the uncertainty on scales far beyond the truncation scales that needs to be considered (e.g., Ricciardulli and Sardeshmukh 2002). Although we have used the Mori–Zwanzig formalism to show that intrinsically there is a need to include the memory term in stochastic parameterizations, there are other ways in the geophysical literature to account for this term. One way basically increases the dimensionality of the system by treating the memory as a red noise (e.g., Newman et al. 1997). Another way is to couple the resolved and unresolved processes explicitly using multiplicative noise. We have used both of these in our study. The memory on resolved time scales has been taken into account by treating it as red noise using the AR(1) process, and the memory on unresolved time scales has been parameterized by the CA-SGS.

Much work remains. We plan to perform a coarse-graining exercise with the FV3GFS model at very high resolution to see if we can optimize the decorrelation length scales chosen for the dynamical memory in our model system. We are also interested in finding out the minimum complexity needed for the dynamical memory, though at the moment we use the first-order autoregressive process AR(1). In the future, we aim to explore the use of cellular automata to represent a wide range of scales, from high-resolution subgrid patterns applied to physical processes to resolutions coarser than the NWP model patterns applied to the tendencies, in order to understand if we can turn off SPPT as our proxy for dynamical memory.

In terms of the physical process perturbations, cellular automata offer a novel approach due to their self-organizing nature. In our application to deep cumulus convection, we found that the model variability of precipitation is improved and that some added variability in the tropics, as well as subtropics, can be seen in the state variables, although there is considerable room for improvement away from the tropics. The impact of the cellular automaton is (not surprisingly) strongly linked to the large-scale state on which we choose to condition it. The choice of large-scale conditioning state gives some room to control the scales, range, and physical meaning of the cellular automata. The framework is quite general, as we can choose which model large-scale state we would like to condition the cellular automaton on and the resolution of the finer grid, depending on the process that we would like to address. Future development includes addressing uncertainty due to turbulence, shallow convection, radiation, and in-cloud microphysics using this general framework.

The stochastically generated skewed (SGS) differential equation is an attractive option to generate physically based perturbations to the field on which we condition the cellular automaton. However, its full advantage for our particular application needs to be further explored and understood. The shape of the distribution given by SGS is linked to the moments given from observations, or large-eddy simulations, and there are challenges deriving a single set of these moments that are applicable to simulating the uncertainty of a subgrid process all over the globe.

## Acknowledgments

This research was partially supported by the Physical Science Division of NOAA’s Earth System Research Laboratory and by the National Weather Service (NWS) Next Generation Global Prediction System (NGGPS) project. The authors would like to acknowledge two NOAA/ESRL/PSD internal reviews by Evelyn Grell and Stefan Tulich, which greatly improved the manuscript. We thank Prashant Sardeshmukh (NOAA ESRL) for helpful discussions on stochastic generated skewed distributions and Rusty Benson (NOAA GFDL) for technical support in the cellular automaton implementation in FV3. We also wish to thank two anonymous reviewers for insightful comments and suggestions leading to further improvements of the manuscript.

## REFERENCES

*Cellular Automata Modeling of Physical Systems*. Cambridge University Press, 341 pp., https://doi.org/10.1017/CBO9780511549755.

*Stochastic Tools in Mathematics and Science*. 3rd ed. Springer, 200 pp.

**72**, 3290, https://doi.org/10.1175/JAS-D-14-0377.1.

*Statistical Mechanics of Nonequilibrium Liquids*. Cambridge University Press, 302 pp.

*Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences*. Springer, 423 pp.

*Nonlinear and Stochastic Climate Dynamics*, C. L. E. Franzke and T. J. O’Kane, Eds., Cambridge University Press, 209–240, https://doi.org/10.1017/9781316339251.009.

*Thermodynamics, Kinetics, and Microphysics of Clouds*. Cambridge University Press, 782 pp.

*Parameterization of Atmospheric Convection. Volume 2: Current Issues and New Theories*, R. S. Plant and J.-I. Yano, Eds., World Scientific, 135–172.

*The Representation of Cumulus Convection in Numerical Models*,

*Meteor. Monogr*., No. 46, Amer. Meteor. Soc., 137–144, https://doi.org/10.1007/978-1-935704-13-3_11.

*Stochastic Numerical Methods: An Introduction for Students and Scientists*. John Wiley & Sons, 381 pp.

## Footnotes

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{1}

A general Langevin equation is a stochastic integro-differential equation relating the particle’s acceleration to its velocity history and to fast fluctuations in the fluid environment (see, e.g., Hohenegger and McKinley 2017). The original Langevin equation describes Brownian motion, which is a special case of the general Langevin equation.

^{2}

Evolution equations including birth and death processes are commonly used to describe the evolution of population characteristics, including population size distributions. Examples for convective cloud population are presented by Plant (2012) and Hagos et al. (2018). Another prominent example is the kinetic equation describing cloud particle size evolution [see, e.g., chapters 13 and 14 in Khvorostyanov and Curry (2014)]. There is no general solution to this equation, but special solutions can be expressed as skewed distributions. In fact, it can be shown that a special solution to the CA evolution equation used in the CA-SGS scheme is a Poisson distribution [see, e.g., chapter 8 in Toral and Colet (2017); Schulman and Seiden (1978)].