Full text
Because the processes controlling Earth's weather and its climatology are complex and inter-related, numerical models are a critical tool for predicting future conditions. Global coverage is necessary because local behavior propagates rapidly to distant areas of the globe. Simulating the whole planet imposes severe computational challenges, however. In the past, this has typically been handled by coarsening model grid spacing until simulations became affordable on the machines of the time. As of 2020, this translated to horizontal grid spacing of ∼10 km for weather models (which simulate days to weeks at a time) and ∼100 km for climate models (which are typically run for centuries). These grid spacings are too coarse to capture many important atmospheric processes.
The impacts of sub-grid scale processes on model climate are instead parameterized based on available grid-scale quantities. Typical parameterized processes include turbulent transport and mixing, gravity-wave motions, greenhouse gas and aerosol chemistry and physics, radiative transfer, and cloud physics. Cloud parameterizations are in particular complicated yet important for accurate predictions. Vapor transport, collisions, and other physics involving micron-scale water drops or ice crystals (collectively called microphysics) are critical for predicting precipitation and future changes in cloud shading. Condensation and evaporation of clouds and resulting fractional cloudiness within a grid cell (often called macrophysics) involve larger spatial scales but are still important to parameterize in conventional models. Condensational heating in convective clouds causes narrow but intense upward vertical motions which are a primary source of vertical transport of heat, moisture, and momentum in the tropical atmosphere (Riehl & Malkus, 1958). Because the microphysics and macrophysics of these intense updrafts are tightly entwined with their motions, convective parameterizations tend to include their own microphysics and macrophysics treatments. Inconsistency between microphysical treatments for convective versus resolved-scale motions is a large source of model biases (Song & Zhang, 2011; Storer et al., 2015). Convection in general has proven to be particularly difficult to parameterize from quantities available on the grid scale (Randall et al., 2003; Stevens & Bony, 2013) and has been implicated as a primary source of climate change uncertainty (Sanderson et al., 2008; Sherwood et al., 2014).
Another challenge posed by coarse resolution is interaction with Earth's surface. Topography is not resolved at typical global model grid spacing and in fact must be even further smoothed to avoid model instability (Lauritzen et al., 2015). Because topography can force air upwards until it condenses, smoothing out high mountain peaks causes major problems for cloud and precipitation climatology (Giorgi & Marinucci, 1996). Insufficient surface roughness means wind stresses are also too weak over smoothed topography and must be parameterized. Subgrid-scale surface heterogeneity also poses problems for coarse models (Prein et al., 2015). And while the focus of this paper is on simulations with prescribed sea surface temperature, it is worth noting that ocean eddies on spatial scales 10 km play a critical role in heat transport (Maslowski et al., 2008) and their parameterization has proven as problematic for ocean models as convective clouds are for atmosphere models (Hewitt et al., 2020). Ocean/atmosphere interaction at convection- and ocean-eddy resolving scales has not (to our knowledge) been studied but is also likely to have important impacts on model behavior.
Because so much is lost at coarse resolution, the global atmospheric modeling community has long pushed toward higher resolution. Unsurprisingly, better topographic resolution improves orographic precipitation, snowpack, and stream flow (Caldwell et al., 2019; Delworth et al., 2012; Duffy et al., 2003; Pope & Stratton, 2002). Sea breeze effects become better captured as coastal boundaries are better resolved (Boyle & Klein, 2010; Love et al., 2011). Because finer grid spacing allows smaller spatial and temporal scales to be resolved, higher-resolution GCMs also better capture extreme precipitation events (Iorio et al., 2004; Terai et al., 2018; Wehner et al., 2014). As GCM grid spacing falls to 25 km or less, tropical cyclones begin to be resolved (Atlas et al., 2005; Bacmeister et al., 2014; Caldwell et al., 2019; Wehner et al., 2014), though capturing details of spatial structure requires still finer resolution (Judt et al., 2021). Some classic model problems are, however, relatively unaffected by reducing grid spacing to 25 km. In particular, increased resolution does not get rid of the erroneous southern branch of the Intertropical Convergence Zone (ITCZ) common in climate models (Bacmeister et al., 2014; Caldwell et al., 2019; McClean et al., 2011). Simulation of the Madden-Julian Oscillation (MJO) is likewise unaffected (Bacmeister et al., 2014; Jung et al., 2012). In addition, precipitation improvement has been found primarily in wintertime (Duffy et al., 2003).
It is notable that these remaining deficiencies are related to convective motions which are unresolved even at high GCM resolutions. Given the aforementioned difficulty of parameterizing convection, this situation is perhaps expected. A small number of global models with grid spacing fine enough to explicitly resolve the largest convection events (hereafter called global convection-permitting models or GCPMs) have also been built. The number of these models has exploded recently because recent advances in computing have tended toward allowing more calculations to be performed in parallel rather than making individual calculations faster. Conventional global simulations already exploit all available parallelism, so won't run faster on these new machines. Higher horizontal resolution is a ready source of increased parallelism, so is attractive in this new computing environment. Unfortunately, smaller timesteps are needed to resolve finer spatial scales. Thus even if all columns could be computed in parallel, a given integration at finer resolution requires more timesteps and therefore has a longer time-to-solution. As a result, GCPM simulations can't be run as routinely nor as long as conventional global models.
The history of GCPM modeling is nicely summarized in Satoh et al. (2019). Briefly, the first GCPM was NICAM, described in Tomita et al. (2005) and Satoh et al. (2008, 2014). For several years its only companion was the Multiscale Modeling Framework (MMF) described in Grabowski and Smolarkiewicz (1999), Randall et al. (2003), and Grabowski (2016). The MMF isn't exactly a GCPM, however, as it replaces the physical parameterizations inside each grid cell of a conventional GCM with a limited-area convection-permitting model (CPM). The MMF is much cheaper than a GCPM because embedded CPMs are typically contained within a single computational node, avoiding MPI communication costs. Additionally, the grid of the CPM is decoupled from that of the GCM, so CPMs are typically 2D and have domain size smaller than the GCM grid cell width. The second GCPM was NASA's GEOS model (Putman & Suarez, 2011), which was used as a synthetic laboratory for designing and testing satellite campaigns (Gelaro et al., 2015) in addition to more general analysis. In the last few years, enough new GCPMs have been developed to warrant their own intercomparison. Called DYnamics of the Atmospheric general circulation Modeled On Nonhydrostatic Domains (DYAMOND), the first phase of this intercomparison focused on a 40-day simulation starting August 1, 2016 and included 8 models with grid spacing less than 5 km globally. An overview of this intercomparison is presented in Stevens et al. (2019). Stevens' study shows striking agreement in outgoing longwave radiation, precipitation, and precipitable water between participating models. Shortwave radiation differs between models, presumably due to differences in low clouds, which aren't well resolved at GCPM resolutions. Models also tend to predict a spurious peak in precipitation just south of the equator, suggesting that km-scale resolution is not the solution to the double-ITCZ problem endemic to conventional climate models (Li & Xie, 2014). Based on the success of this first intercomparison, a second DYAMOND intercomparison (called DYAMOND2) is now underway. The current paper documents a new contribution to DYAMOND2.
GCPMs can be viewed as a natural and beneficial extension of conventional GCMs to finer resolution, but they can also be seen as the extension toward larger domains of a robust research community focused on limited-area CPMs. Beginning with the explicit simulation of a single convective event (Ogura, 1963), cloud-resolving simulations have steadily grown in duration and domain size. Recently, Bretherton and Khairoutdinov (2015) and Narenpitak et al. (2017) describe multi-month 4 km simulations simulating the entire tropical channel between 45°N and 45°S. CPMs tend to offer more benefit for summertime convection rather than wintertime cyclones (Prein et al., 2015), as may be expected given the spatial scale of these storm types. Limited-area CPM research suggests that resolution finer than ∼4 km is needed to resolve convective ensemble statistics (Weisman et al., 1997; also found for GCPMs by Miyamoto et al. [2013]) but resolution finer than that adds relatively little value (Kain et al., 2008; Langhans et al., 2013; Schwartz et al., 2009). Cloud fraction tends to decrease as resolution becomes finer (Fosser et al., 2014; Langhans et al., 2013; Prein et al., 2013), a feature also found in GCPMs (Hohenegger et al., 2020; Noda et al., 2010).
A great deal of CPM research has been organized around the Global Energy and Water Cycle Experiment Cloud Systems Study (GCSS). As described in a review by Krueger et al. (2016), GCSS organized intercomparisons of CPMs and single-column versions of GCMs for intensive observing periods spanning a wide variety of cloud regimes. These intercomparisons clarified processes CPMs could and couldn't handle, often leading to idealized follow-up experiments. These follow-up studies have proven invaluable for providing process insights and subsequent model improvements. DYAMOND is in some ways the reincarnation of GCSS for the next generation of models.
In general, high-resolution regional studies have added value primarily by resolving fine-scale features rather than through upscale effects onto scales resolved by conventional models (Caldwell, 2010; Prein et al., 2015). One potential reason for this is that lateral boundary conditions impose strong constraints on domain-averaged properties (Edman & Romps, 2014). Thus while GCPMs may be overkill for looking at fine-scale features which could be studied via limited-area models, they offer fresh new potential to solve long-standing deficiencies in the general circulation.
The goal of this paper is to introduce the GCPM being developed by the Energy Exascale Earth System Model (E3SM) project and to provide an initial look at its behavior in the DYAMOND2 case study. Details about this model are provided in Section 2. Sections 3–5 describe experimental design, data for evaluation, and computational performance (respectively). Results in Section 6 are broken into an analysis of effective resolution in Section 6.1, general attributes in Section 6.2, clouds and radiation in Section 6.3, precipitation in Section 6.4, and specific weather phenomena in subsequent subsections. Conclusions follow in Section 7.
Model DescriptionAs described in Golaz et al. (2019), the E3SM project was born from the US Department of Energy (DOE)'s need for quantitative information about future climate for use in energy-sector decisions. Given DOE's leadership in high-performance computing, it has been natural for E3SM to focus on compute-intensive frontiers in climate science. One of those efforts has been to develop a new GCPM called the Simple Cloud-Resolving E3SM Atmosphere Model (SCREAM).
Our ultimate goal is to make SCREAM as fast as possible on exascale machines by writing it in C++ using the Kokkos library (Carter-Edwards et al., 2014) for performance portability. See Bertagna et al. (2019) and Bertagna et al. (2020) for a description of our design strategy and initial performance results. We are, however, approaching this goal by first creating a prototype version in Fortran using the existing E3SM atmosphere infrastructure. This initial implementation—which is the focus of the current study—is being used as the template for the C++ implementation as well as giving us an early look at model behavior. The final implementation should be scientifically identical to this prototype version but will be much faster because of its ability to run on GPU-powered computers.
Our strategy has been to make our first implementation as simple as possible and to start using it for science as quickly as possible. This strategy is expected to result in sub-optimal skill in our first implementation, but allows us to more rapidly produce, understand, and improve our model. We believe that it is better to start with an overly simple model and to add complexity as needed rather than to start with a more sophisticated/accurate model which we don't understand.
Simplicity in particular means that SCREAM consists solely of nonhydrostatic fluid dynamics, a turbulence/cloud fraction scheme, a microphysics scheme, a radiation scheme, an energy fixer, and prescribed-aerosol functionality. These pieces are described in the subsections below. SCREAM does not parameterize sub-grid scale gravity-wave drag or deep convection. This initial implementation uses the E3SM land model described in Golaz et al. (2019). It also uses prescribed-ice mode from CICE4 (Hunke & Lipscomb, 2008) to compute surface fluxes, snow depth, albedos, and surface temperature, resetting sea ice thickness after each timestep to 2 m in the northern hemisphere and 1 m in the southern hemisphere. Sea surface temperature (SST) is prescribed.
Fluid DynamicsSCREAM's fluid-dynamics solver (hereafter dycore) solves the nonhydrostatic equations of motion in a rotating reference frame with the shallow atmosphere approximation and a hyperviscosity based turbulence closure. It additionally transports several constituents, including multiple forms of water and various aerosols. It is implemented in the High Order Method Modeling Environment (HOMME) (Dennis et al., 2005, 2012; Evans et al., 2013). HOMME contains several spectral element based dycores, including the hydrostatic dycore used by E3SM (Caldwell et al., 2019; Golaz et al., 2019; Rasch et al., 2019) and the Community Earth System Model (Small et al., 2014; S. Zhang et al., 2020). We refer to the new nonhydrostatic dycore developed for SCREAM as HOMME-NH.
HOMME-NH uses the nonhydrostatic formulation of the equations from Taylor et al. (2020). It solves the equations in a terrain following mass based vertical coordinate (Kasahara, 1974; Laprise, 1992), with prognostic equations for the three components of the velocity field, the mass-coordinate pseudo-density, the geopotential height, and a thermodynamic variable, for which we use virtual potential temperature. The prognostic equations consist of the time-reversible adiabatic terms from Taylor et al. (2020), a ∇4 hyperviscosity following Dennis et al. (2012) and Guba et al. (2014), and a sponge layer at the model top. For the adiabatic terms, we use a structure preserving formulation in order to preserve the discrete Hamiltonian and produce an energetically consistent model. The horizontal discretization uses the collocated mimetic spectral element method from Taylor and Fournier (2010), with conservative and monotone semi-Lagrangian tracer transport (Bradley et al., 2019). The vertical discretization uses a Lorenz staggered extension of the mimetic centered difference from Simmons and Burridge (1981). With this vertical staggering, prognostic variables are located at level midpoints, with the exception of the vertical velocity and geopotential, which are located at level interfaces. For the vertical transport terms, we use a vertically Lagrangian approach adapted from Lin (2004).
For the temporal discretization, we use a Horizontally Explicit Vertically Implicit (HEVI) approach (Satoh, 2002), discretized with an IMplicit-EXplicit (IMEX) Runge Kutta method (Ascher et al., 1997). The HEVI splitting decomposes the equations into a set of terms which represent vertically propagating acoustic waves (treated implicitly), and the remaining terms which include all horizontal derivatives (treated explicitly). We use a highly efficient IMEX method from Steyer et al. (2019) and Guba et al. (2020), with a second-order accurate coupling of a high-stage high-CFL scheme for the explicit terms and a Diagonally Implicit Runge Kutta (DIRK) scheme for the implicit terms. Due to the use of the Laprise mass coordinate, the vertical acoustic waves are isolated to only two terms in the equations for vertical velocity and geopotential solved at level interfaces, leading to an implicit system for a single variable.
There are several sources of dissipation in the dynamical core. The ∇4 hyperviscosity is the largest. It is applied to all prognostic variables and on every model layer, with a hyperviscosity coefficient of 2.5 × 1010 m4 s−1 for the 3.25 km grid. Because tuning at 3.25 km is expensive, we chose this value based on a Δx3 scaling of the hyperviscosity coefficient used by E3SM at lower resolutions. For the model-top sponge-layer, we applied a ∇2 Laplacian operator to all prognostic variables according to the reference-pressure based ramp function from Lauritzen et al. (2011). This ramp starts at layer 14 (∼19 hPa) with a coefficient of 0.189 × 10−4 m2 s−1 and ramps up to 6.93 × 10−4 m2 s−1 at the model top. In addition, vertical dissipation is introduced by the monotone vertical remap operator. A smaller amount of dissipation is also generated by the Runge-Kutta timestepping.
Model GridOur horizontal grid for dynamics is a cubed-sphere grid with 1,024 × 1,024 spectral elements on each face, denoted ne1024. The total number of elements is therefore 6,291,456. Within each element, fields are represented by degree-3 polynomials, using nodal values on a 4 × 4 grid of Gauss-Lobatto-Legendre (GLL) nodes. The edge and corner nodes are shared by adjacent elements, resulting in an average spacing between GLL nodes of ∼3.25 km. The nonuniform spacing of GLL nodes presents some challenges to the physical parameterizations (Herrington et al., 2019), which we avoid by evaluating the parameterizations on a uniformly spaced 2 × 2 grid within each spectral element. This physics grid has 4/9 as many physics columns as would be in a GLL-collocated physics grid. Tests show that the 2 × 2 physics grid provides very similar results to simulations with physics running on every GLL node (Hannah et al., 2021). Our land model is run on a 1/8° latitude-longitude grid. SST and sea ice are computed on the high-resolution ocean grid used by Caldwell et al. (2019), which tapers from 18 km in the tropics to 6 km near the poles. The ocean and sea ice grids have minor impact since SST and ice extent are interpolated from 0.5° data sets. It would be better to have all surface calculation on the 3.25 km atmosphere grid, but resolution challenges with the E3SM input data toolchain made doing so impractical for this initial simulation.
We use a relatively fine 128 layer vertical grid with a model top at 40 km (2.25 hPa) and a sponge layer above ∼19 hPa (as as described in the previous subsection). Vertical grid spacing is presented in Figure 1. Representative grid spacing in the boundary layer is ∼50 m, in trade Cu is ∼100 m, and in tropical cirrus anvils is ∼250 m.
Enlarge this image.Figure 1. SCREAM grid spacing. Panel (a) shows the complete vertical grid using logarithmic pressure to emphasize the upper atmosphere. Panel (b) zooms in on the troposphere using linear pressure spacing to emphasize lower levels. The sponge layer is indicated by gray-shading.
To generate the SCREAM v0 surface topography, we use the NCAR topography tool chain (Lauritzen et al., 2015) to first compute the unfiltered height field on the atmosphere grid. We then smooth the height field on the GLL grid using 16 iterations of the spectral element Laplace operator. To quantify the amount of smoothing, we follow Evans et al. (2013) and compare power spectra E(k) from the spherical harmonic representation of the filtered and unfiltered height field, and then compute the lowest wave number k1/2 for which the smoothing has reduced E(k1/2) by 50%. The SCREAM v0 topography has a k1/2 corresponding to wavelength 6.4Δx.
Clouds and TurbulenceBoundary layer clouds and their associated circulations are still largely unresolved at 3.25 km so a parameterization of interaction between clouds and turbulence is critical. Because GCPMs push the boundary of computational possibility, it is important that these processes are handled efficiently. These goals are accomplished in SCREAM via the Simplified Higher Order Closure (SHOC; Bogenschutz & Krueger, 2013). Similar to other widely used assumed PDF-based schemes (Cheng & Xu, 2008; Golaz et al., 2002), SHOC computes subgrid-scale liquid cloud and turbulence using an assumed double-Gaussian probability density function (PDF). SHOC is more efficient than the aforementioned schemes, however, because it diagnoses rather than prognoses the higher order moments that are needed to close the double Gaussian PDF. Bogenschutz and Krueger (2013) demonstrate that when SHOC is used in limited-area cloud-resolving simulations of boundary layer clouds, the solution is insensitive to the horizontal resolution choice. This is in contrast to a standard 1.5-order TKE closure, which suffers from large horizontal-resolution sensitivity when used in the same cloud-resolving model.
SHOC has undergone several updates since Bogenschutz and Krueger (2013) to improve numerical stability and performance among the wider range of regimes SHOC is subjected to in a global model. Chief among these updates is the implementation of an implicit diffusion solver, a revised formulation of the turbulence length scale to better achieve vertical convergence, and a revised formulation of the eddy diffusivities for the stable boundary layer (similar to those implemented in Bretherton and Park [2009]). The turbulence length scale is now a continuous formulation that avoids the separate definitions of in-cloud versus sub-cloud length scales documented in Bogenschutz and Krueger (2013) and performs scientifically similarly to the original formulation.
In addition to the liquid cloud fraction supplied by SHOC, we require an ice cloud fraction. For simplicity, our initial implementation includes the same ice cloud fraction used by E3SMv1 and inherited from CESM1. This implementation assumes ice cloud starts forming when an ice-modified relative humidity RHi = (qv + qi)/qsat, i reaches a user-specified minimum value and reaches 100% at a user-specified maximum value. Unfortunately, these parameters were left at their low-resolution E3SMv1 defaults of 80% and 105% (respectively) in our DYAMOND2 simulation. The impact of this mistake is shown in Section 6.2.
MicrophysicsSCREAM microphysics is based on the Predicted Particle Properties (P3) scheme of Morrison and Milbrandt (2015) taken from version 4.1 of the Weather Research and Forecasting (WRF) model (Skamarock et al., 2019). The novel feature of P3 is that it avoids arbitrary cutoffs between cloud-borne and precipitating ice categories by employing a single ice category which is allowed to evolve naturally from small pristine crystals into large and possibly rimed snowflakes. While the WRF version of P3 allows for multiple simultaneous populations of these ice crystals within a grid cell, SCREAM currently only supports a single population because the modest improvements from multiple ice populations reported in Milbrandt and Morrison (2016) were not deemed worth the additional software engineering time required to support this feature. The liquid phase of the P3 scheme—like most microphysics codes—separates drops into cloud and rain categories.
One feature of this scheme is the clever use of supersaturation to diagnose condensation, evaporation, sublimation, and deposition. This approach works well for Large-Eddy simulations (LES) which explicitly model each updraft, but probably underpredicts condensation for the 3.25 km grid spacing used in SCREAM (Morrison & Grabowski, 2008). The great benefit of this supersaturation approach is that it treats ice growth at the expense of nearby liquid (Bergeron, 1935; Findeisen, 1938; Wegener, 1911, hereafter WBF process) in a very natural way. Unfortunately, allowing supersaturation in P3 directly conflicts with the instantaneous saturation adjustment assumption which forms the foundation of SHOC's PDF. For consistency, our P3 implementation instead handles vapor deposition, sublimation, and the associated WBF process following Gettelman and Morrison (2015). In particular, maximum overlap between liquid and ice is assumed when liquid and ice coexist, leading to efficient liquid-to-ice transition via the WBF process. If all liquid is removed within a microphysics timestep, vapor deposition onto ice for the remainder of that timestep is computed based on cell-average water vapor content.
Another inconsistency between SHOC and the WRF version of P3 is the use of fractional cloudiness and precipitation. P3 neglected all sub-grid variability such that cloud and precipitation covered the entire grid cell where they exist and otherwise the cell was entirely devoid of condensate. SHOC provides fractional cloudiness, so we modified P3 to only operate in the cloudy or precipitating portion of each cell. Our fractional cloudiness implementation is similar to Jouan et al. (2020), which was implemented in WRF P3 around the same time as we made our modifications. The fraction of each cell containing precipitation is also important. In SCREAM this was taken to be equal to the largest cloud fraction of all cells including and above the layer of interest. This approach is crude (as noted by Zheng et al. [2020]) and will be a subject of future research.
SHOC's subgrid assumptions require further modifications. SHOC uses a double-Gaussian PDF to model subgrid-scale variations in liquid water potential temperature, total water mixing ratio, and vertical velocity. Larson and Griffin (2013) provide an analytical formulation for incorporating SHOC's variability into microphysical processes expressed as power functions. We intend to implement this consistent scheme in our version of P3 eventually, but for the moment we have instead implemented the partially consistent approach from Morrison and Gettelman (2008), which instead assumes a gamma distribution for liquid water mixing ratio and ignores subgrid temperature variations. The benefit of the gamma distribution is that the expected value of a power-law-based microphysical process rate can be written as that power law applied to the cell-mean value multiplied by an easily calculated scaling factor.
Finally, water vapor saturation was changed in our version of P3 to be consistent with the Murphy and Koop (2005) (MK) implementation used in SHOC. MK is more accurate at very low temperatures than the Flatau et al. (1992) implementation originally used in P3, but is more computationally expensive. We found this performance difference, however, to have a negligible impact on total run time.
RadiationGas optical properties and radiative fluxes are computed using the RTE + RRTMGP radiative transfer package (Pincus et al., 2019). Active gases in SCREAM include H2O, CO2, O3, N2O, CO, CH4, O2, and N2. Cloud and aerosol optical properties are computed as in the Community Atmosphere Model (CAM). The approach is described in detail in Neale et al. (2012). Briefly, condensed phase optical properties (extinction coefficient, single scattering albedo, and asymmetry parameter for shortwave bands and absorption coefficient for longwave bands) are computed per unit mass for liquid, ice, and aerosol, then multiplied by the appropriate mass mixing ratio for use in RTE + RRTMGP.
Liquid cloud optical properties are calculated from a table-lookup after being computed offline using a Mie scattering code (Wiscombe, 1996) based on the assumption (taken from microphysics) that the total number of liquid drops with diameter D follows a gamma histogram [Image Omitted. See PDF]
with intercept parameter N0, slope parameter λ, and spectral size dispersion μ taken every timestep from P3. In this initial implementation, in-cloud liquid water content is assumed to be homogeneously distributed. This is inconsistent with our implementation of P3, which (as noted above) assumes a gamma distribution for spatial variations in cloud liquid. Fixing this inconsistency is a future goal.
Ice cloud optical properties are computed for each shortwave and longwave band used by the radiation code using a lookup table based on the modified anomalous diffraction approximation (Mitchell, 2002). The only input to these table lookups is ice effective radius, which is computed in P3. Because ice mass-density relationships are different for different size and riming regimes, ice effective radius is calculated via a table lookup described in Morrison and Milbrandt (2015). Because P3 merges the ice and snow categories used by traditional microphysics schemes into a single ice mode, radiation naturally acts on all frozen hydrometeors. Aerosol optical properties are specified in a lookup table as a function of wet refractive index and wet surface mode radius (Ghan & Zaveri, 2007).
Vertical overlap of partially cloudy cells is accounted for by assuming maximum-random overlap (Geleyn & Hollingsworth, 1979) using the Monte Carlo Independent Column approach (MCICA; Pincus et al., 2003).
Prescribed AerosolE3SMv1 uses a 4 Mode Aerosol Model (MAM4 Liu et al., 2016). For computational efficiency, we employ a version where this modal aerosol information is prescribed using monthly average climatologies interpolated to the model grid from a 1° resolution E3SMv1 simulation. Implementation and use of prescribed-aerosol functionality is described in K. Zhang et al. (2013), Lebassi-Habtezion and Caldwell (2015), and Shi and Liu (2018). The default prescribed-aerosol implementation scales aerosols by different random perturbations every day to improve agreement between prescribed- and prognostic-aerosol simulations at high latitudes. These random daily jumps are confusing for analysis of short timeseries, so we've set the magnitude of random perturbations to zero for DYAMOND2. This might degrade aerosol behavior in polar regions.
Like E3SMv1, cloud condensation nuclei (CCN) concentration is derived from Abdul-Razzak and Ghan (2000). Ice nucleation follows Gettelman et al. (2010) for deposition nucleation and homogeneous freezing of solution droplets but retains the original P3 implementation for cloud and rain drop freezing.
Energy FixerSCREAM inherited its energy fixer from CAM. As described in Lauritzen and Williamson (2019), this energy fixer corrects errors due to pressure work, time integration in the dynamical core, inconsistent formulations of equation of state, and other minor sources of non-conservation. Historically, CAM and the atmospheric component of E3SM had used an incorrect formulation for energy. Williamson et al. (2015) documents this problem and provides a correction, which is used in SCREAM.
TimestepsLike most atmosphere models, SCREAM's timestepping is a complex mixture of substepping and superstepping of individual processes. Ideally, model timesteps would be small enough that modest changes wouldn't have a noticeable effect on model behavior. Unfortunately, climate models have not yet reached that goal (Santos et al., 2020). Thus we list the timesteps used for the DYAMOND2 simulation in Table 1.
Table 1 Timesteps Used in SCREAM DYAMOND2 Simulation (in sec)
| Main | Dycore | Dycore Remap | Advection | Radiation |
| 75 | 9.375 | 18.75 | 75 | 300 |
Note. Processes not listed use Main timestep.
TuningTuning is important for optimal performance of any weather or climate model, but should become less important at higher resolution where more processes are explicitly resolved and therefore expressed in a more complete and physical way. Because of time constraints and a reticence to tune away problems before understanding their source, the only parameter adjustment we made was to modify the lower limit of the eddy diffusivity damping timescale to get net top-of-atmosphere (TOA) radiation to match observations and to control surface temperatures under stable conditions at high latitudes. Because our tuning was based on short (1 or 2 days) simulations and therefore required comparison against higher-time frequency radiative observations which (as described in Section 4) have larger global-average bias than the monthly average data used to assess the simulation-average radiation, the TOA net bias reported here still ended up being somewhat large. Our crude tuning approach also resulted in clouds which are too stratiform rather than convective (as described in Section 6.3). High latitude land surface temperature biases remain high, indicating that more tuning work is needed.
Experimental DesignThe focus of this study is a 40-day global simulation (January 20 through February 28) performed as part of the DYAMOND2 intercomparison. Our implementation follows the guidance at
As mentioned in Section 2.7, aerosol distributions are prescribed from a 1° E3SMv1 simulation. This simulation was 6 years long with annually repeating forcings (SST, sea ice extent, land use, solar forcing, aerosol emissions, greenhouse gases, and volcanic aerosols) values typical for the decade surrounding 2010. The last 5 years of this simulation are averaged to create a monthly varying aerosol field.
Soil and snowpack initial conditions were computed in 2 steps. First, the E3SM land model was run from January 1, 1979 through August 1, 2016 at the target resolution forced by observed atmospheric conditions from Version 7 of the Climatic Research Unit—National Centers for Environmental Prediction (CRUNCEPv7, Viovy [2018]) atmospheric forcing data. This simulation couldn't be extended beyond 2016 because of CRUNCEPv7 data availability. The second step was therefore to run from August 1, 2016 to January 20, 2020 using EAMv1 at 1° nudged to ERA5 reanalysis with a 6 hr timescale. Prescribed weekly SST and sea ice from OISSTv2 (Reynolds et al., 2002) is used for this simulation. The machinery for this second step came from the Cloud-Associated Parameterizations Testbed (Ma et al., 2015; Phillips et al., 2004).
Observations for EvaluationThe short duration of this simulation and our focus on small time and spatial scales limit the range of observational data sets suitable for comparison. We rely heavily on the European Center for Medium-Range Weather Forecasting's ERA5 reanalysis (Hersbach et al., 2020). This retrospective simulation assimilates a massive array of observations, runs at 31 km horizontal resolution with 137 vertical levels and a top at 0.01 hPa, and is available at hourly resolution. Because model formulation strongly affects cloud and precipitation predictions from reanalysis, we use satellite products for cloud-related variables. In particular, we use half-hourly 0.1° gridded Global Precipitation Measurement (Hou et al., 2014, GPM) Integrated Multi-satellitE Retrievals for GPM (IMERG) product version V06 B (G. J. Huffman et al., 2019) for global precipitation. For radiative fluxes, we use CERES-EBAF 1° data averaged over February 2020 (Loeb et al., 2018). To examine the radiative properties of individual storms, we also use CERES-SYN hourly 1° data (Doelling et al., 2013, 2016). Cloud fraction and liquid water content are taken from CloudSat (Austin et al., 2009; Su et al., 2011) and from the CERES–CALIPSO–CloudSat–MODIS merged product (Kato et al., 2010, C3M). CloudSat and C3M are not available for the 2020 dates simulated and are instead climatological averages.
Where possible, we compute long-term averages using the last 30 days of the simulation (January 30 through February 28); we exclude the first 10 days of the run as spinup (though SCREAM fields stabilize after just one day of spinup, see Figure 5). As noted above, some observational data sets are only available as monthly averages. For corresponding variables, we show results using just days in Feb. Finally, the first week or so of the simulation can be treated as a weather forecast, we use all 40 days of the simulation for some analysis of storm behavior.
PerformanceThe DYAMOND2 simulation was performed as a series of 1536-node job submissions using the Knights Landing (KNL) nodes of Cori at the National Energy Research Supercomputing Center (NERSC). We found that using 8 MPI processes and 16 OpenMP threads per node provided the optimal balance of memory usage and performance for these 1536-node jobs. The overall throughput for the 40-day simulation, including I/O, was about 4–5 simulated days per day (SDPD). Further details about the performance of this 40-day DYAMOND2 simulation can be explored at
Table 2 SCREAM Timings as a Function of KNL Node Count Using Either 8 × 16 MPI Tasks Versus OpenMP Threads or 16 × 8 MPI Tasks Versus OpenMP Threads Per Node
| Nodes | 8 × 16 Dycore timing in minutes | 8 × 16 with IO SDPD | 8 × 16 without IO SDPD | 16 × 8 without IO SDPD |
| 1536 | 100.8 | 5.1 | 5.8 | OOM |
| 3072 | 53.9 | 8.6 | 10.3 | OOM |
| 4096 | 44.4 | Not run | Not run | 14.2 |
| 6144 | 29.2 | 14.2 | 19.2 | 23.1 |
Note. All timing runs were 1 day in length. Timings with IO include all standard output for our DYAMOND simulation. OOM means Out of Memory and IO stands for Input/Output.
The simulation used the Software for Caching Output and Reads for Parallel I/O (SCORPIO) library for reading input data and writing simulation output to the file system. SCORPIO is derived from the Parallel I/O library (Hartnett & Edwards, 2021) and continues to support the same application programming interface. To improve the I/O write performance the library caches and rearranges output data between MPI processes before using low level I/O libraries like the netCDF, Parallel netCDF (PnetCDF) (Latham et al., 2003), and ADIOS (Godoy et al., 2020) libraries to write the data to the file system. On Cori the simulation produced ∼4.5 TB of data per simulated day and achieved an average I/O write throughput of ∼2.5 GB/s using the PnetCDF library.
Unsurprisingly for such a large run, we experienced several node failures during the simulation requiring restarts from the previous day. Because E3SM is bit-for-bit reproducible for identical initial conditions and forcings, these failures should not have any impact on our results. During model development, we had problems with occasional extremely cold temperatures near the surface at wintertime high latitudes. We fixed this problem by increasing turbulent diffusivity in stable atmospheric conditions, but this had the side effect of increasing time-average warm bias in polar regions. The tuning used here balances model stability against bias.
Results Kinetic Energy SpectrumAt convection permitting resolutions, the simulated atmosphere's kinetic energy spectra recovers many features seen in observations and reveals many aspects of model diffusion, filtering and parameterization behavior (Skamarock et al., 2014). As a first look at this in SCREAM, we plot the horizontal kinetic energy power spectra at 250 and 500 hPa in Figure 2. The spectra are computed via spherical harmonic transforms of 3-hr flow snapshots from days 22 to 23 of the simulation. We denote by E(k) the power of the spherical harmonics of degree k. We plot compensated spectra, E(k)k5/3, to better illustrate the high wave number k−5/3 regime. SCREAM reproduces the observed Nastrom-Gage transition from a k−3 scaling at low wavenumbers to a k−5/3 regime (Lindborg, 1999; Nastrom & Gage, 1985). The k−5/3 region extends to ∼6Δx wavelength (wavenumber 2000), where the spectra start to roll-off and become dominated by model diffusion. Thus SCREAM's effective resolution is similar to ICON and IFS (Neumann et al., 2019) despite SCREAM's novel use of a coarser grid for physical parameterizations.
Enlarge this image.Figure 2. Compensated kinetic energy spectra (E(k)k5/3) at 500 and 250 hPa from days 23 to 24 of the simulation. The black lines show idealized E(k) ≈ k−3 and E(k) ≈ k−5/3 scalings. See text for details.
Global-average model biases are modest in size but are generally larger than the range of observed day-to-day variability within the simulation period (Figure 3). TOA net shortwave (SW) radiative absorption SWnet and longwave (LW) emission LWnet are both too strong but (as noted in Section 2.10) were tuned to compensate each other such that TOA radiative bias radnet exhibits only a modest warming tendency. Radiative biases are almost entirely due to clouds rather than clear-sky bias (not shown). Too little SWnet reflection and excessive LWnet emission suggests a lack of clouds, so it is surprising that model calculated vertically projected cloud fraction is 5% too large. This is an unfortunate result of using a RH-based ice cloud fraction parameterization without retuning for higher resolution. As a result, large cloud fraction occurs in cold regions which don't necessarily have cloud mass (Figure 4). Fortunately, “clouds” without condensate are treated like clear-sky air by radiation, so our mistake is mostly cosmetic in nature. In the future we intend to switch to a mass-based all-or-nothing ice cloud fraction scheme to avoid this problem. An offline version of this mass-based approach (shown in Figure 4) is used in the remainder of this paper wherever upper-level cloud fraction is required.
Enlarge this image.Figure 3. Global-mean anomaly in variables listed along x-axis. Anomalies are calculated relative to the February 2020 average of CERES-EBAF (for radiative fluxes and cloud fraction), and the January 30, 2020 through February 28, 2020 average of ERA5 (for precipitable water, sensible and latent heat fluxes, and near surface temperature), and GPM (for precipitation). Each dot represents a single daily average, so vertical spread gives a sense of temporal variability. There is a dot for each day in February 2020. Units for each variable are included in the x-axis labels.
Enlarge this image.Figure 4. Vertical profile of Feb-mean tropics-averaged (30°S–30°N) cloud fraction computed by SCREAM compared to an offline calculation of cloud fraction based on assuming an entire cell is saturated whenever cloud water content >10−5 kg kg−1.
Global-average precipitation is ∼0.3 mm day−1 larger in SCREAM than GPM, which is consistent with a general tendency for models to have higher precipitation rates than observations (Terai et al., 2018), including in the previous DYAMOND intercomparison (Stevens et al., 2019). Temperature at 2 m height (T2m) and vertically integrated vapor lie within observed day-to-day variability in the global average, though we show later that this is due in part to compensating errors. Sensible heat flux (SHF) and surface evaporation (a.k.a. latent heat flux; LHF) are larger than observed, probably due to near-surface wind speed biases discussed later.
Figure 5 demonstrates that our simulation doesn't drift rapidly in time, even in the first few days of the run. Time tendencies in other key variables are likewise small (not shown). The amplitude of global-average diurnal variations is also reasonable. Interestingly, GPM and ERA5 contain periods where global-average precipitation drops, while SCREAM is more temporally invariant. Understanding what causes these global-average drops is an interesting question for future work.
Enlarge this image.Figure 5. Fifteen minute (thin curves) and daily mean (thick curves) time series of global-average precipitation (left) and 2 m temperature (right) for the duration of the DYAMOND2 simulation.
Near-surface temperature biases are modest at low latitudes and larger at high latitudes (Figure 6). In the first few days of our simulation, T2m was uniformly too high at high latitudes (not shown), which we attribute to a land initial condition created by driving our land model with a 1° atmosphere model which one might expect to handle snowpack poorly. We tuned overturning turbulent mixing in stable conditions to compensate the warm biases we saw in our initial short testing runs; it appears in retrospect that we overdid it. Averaged over the last 30 days of the simulation, the US, Greenland, and the far eastern side of Russia retain 6 K warm biases, while north Asia and the Canadian Arctic are ∼5 K too cold. Improving these temperature biases is a future goal. Surface pressure is also too large at high latitudes (right-hand panel of Figure 7), which will translate (through thermal wind balance) to errors in wind speed.
Enlarge this image.Figure 6. Near-surface temperature averaged over January 30 through February 28, 2020 from SCREAM and ERA5 reanalysis.
Enlarge this image.Figure 7. Bias (relative to ERA5) in 10 m wind speed (left) and surface pressure (right) averaged over the last 30 days of the simulation.
Near-surface wind speed is too high almost everywhere but particularly over midlatitude oceans (Figure 7). Bias is smallest in the tropics. We are still working to understand and fix this deficiency, but note that switching to the Zeng et al. (2002) scheme significantly alleviates excessive wind speeds. Consistently positive wind biass is solely a feature of the surface layer—even at 925 hPa wind biases are much more balanced around zero. Overall, it is surprising that so many aspects of our simulation look quite good in spite of this near-surface wind bias. Overly strong SHF and LHF mentioned earlier are unsurprising given strong near-surface wind speed.
Figure 8 shows geopotential height and wind speeds on the 200 hPa pressure surface averaged over the period from January 30 to February 28. Although there is generally strong agreement between SCREAM and ERA5, two hotspots emerge. First, the wintertime Rossby wave train that reinforces the upper-level trough over Greenland is markedly more intense in SCREAM than in ERA5. The result is southward displacement of the subtropical jet (STJ) over the West Atlantic and anomalously strong poleward flow from the STJ toward Greenland. In fact, this anomaly in the Central Atlantic is largely barotropic, present even at 850 hPa with approximately the same magnitude (not shown). A second region of anomalous behavior also exists around the periphery of Australia where the 200 hPa geopotential surface is enhanced, producing spurious meridional flow throughout this region. Notably, the bias pattern present in the difference plots suggest an enhancement in wavenumber 4 in both hemispheres centered around the locations of cubed-sphere corners in the dynamics grid. The bias appears slightly stronger in the first 20 days of the simulation than the last 20 days (not shown). The source of this behavior is under investigation.
Enlarge this image.Figure 8. 200 hPa geopotential height (top), zonal wind speeds (middle), and meridional wind speeds (bottom) averaged over the last 30 days of the simulation. Stippling in the difference plots (right panels) indicates regions where SCREAM falls outside the range of mean values for all years in ERA5 1979–2020.
SWnet and LWnet radiation biases were found in Figure 3 to somewhat cancel in the global mean; Figure 9 reveals that this cancellation also holds regionally in many places. Cancellation between SW and LW biases is a hallmark of high clouds. Further evidence of problems with high clouds is the pattern of LWnet bias, which is large where deep convective clouds are expected.
Enlarge this image.Figure 9. TOA radiation averaged over February 2020. Top is SW (>0 warms the planet), middle is LW (>0 cools the planet), and bottom is net (>0 warms the planet).
Figure 10 explores the vertical profile of tropical clouds compared to climatological CloudSat measurements. Because SCREAM results are for one month only, detailed comparison is not appropriate. Nonetheless, SCREAM's ability to capture the general features from CloudSat data is very good, particularly compared to the (albeit old) GCMs analyzed in Su et al. (2011). In particular, SCREAM captures the bimodality of deep and shallow clouds and does a reasonable job of matching the quantitative magnitude of each peak. Ability to better capture the structure of tropical convection is perhaps unsurprising given that resolving such convection was a primary motivation for developing a 3.25 km model. Both simulated cloud peaks sit lower in the atmosphere than they do in the measurements. Another notable deficiency in SCREAM is the lack of mid-level clouds, which may be tied to either the absence of significant cloud detrainment at mid-levels, overly efficient sedimentation of cloud particles through mid-layers, or both. Reasonable or even excessive SCREAM anvil condensate in Figure 10 and erroneously large high cloud fraction in Figure 4 are at odds with excessive LW emission to space in Figure 9. We are still working to understand this conundrum.
Enlarge this image.Figure 10. Cloud water content (CWC) profiles from SCREAM (solid) versus CloudSat observations from Su et al. (2011) (dashed). Data are averaged over all longitudes and latitudes between 30°S–30°N and over all 40 simulated days.
Net outgoing radiation over the northern hemisphere oceans is found in Figure 9 to be too strong in general (i.e., the oceans in Figure 9i are colored blue indicating more radiation leaving than entering the atmosphere). This is due mainly to trapping of LW emission; SWnet insolation at higher northern latitudes is too small in wintertime to matter.
Away from high-latitude winter regions, the impacts of high clouds on SWnet and LWnet tend to cancel so radnet is a good indicator of lower-level cloudiness. Figure 9 reveals a lack of low clouds over the southern ocean, but generally decent low-cloud radiative forcing in the stratocumulus decks off the west coast of the continents. Anemic stratocumulus is a perennial GCM bias (Jian et al., 2020; Nam et al., 2012), so capturing this cloud type in SCREAM is exciting. This is particularly surprising since 3.25 km grid spacing is generally considered insufficient to capture boundary-layer clouds like this. One potential reason for improvement is our higher-order turbulence closure. Increased vertical resolution (∼50 m in the boundary layer) in addition to SCREAM's high horizontal resolution also likely helps; Bogenschutz et al. (2021) and Lee et al. (2021) demonstrate that increased vertical resolution helps to ameliorate these biases in E3SM, owing to better representation of the cloud top cooling and turbulence feedback, but both studies hypothesize that concurrent increases in the horizontal and vertical resolution are needed to adequately simulate the coastal Sc. Results with SCREAM support that hypothesis.
Figures 11a and 11b display the February 2020 average profiles of cloud fraction and cloud liquid water for SCREAM and the February 2006–2010 climatology from C3M. These profiles are averaged over a small domain neighboring the coast of Peru and Chile. This domain was selected as it represents the area of most intense shortwave cloud radiative effect (SWCRE) biases associated with low clouds in the northern-hemisphere winter season for standard-resolution GCMs (e.g., Danabasoglu et al., 2020; Golaz et al., 2019). Although different averaging periods are used for C3M versus SCREAM data, stratocumulus are a persistent feature in this region so broad comparison is reasonable. SCREAM produces cloud structure quite similar to the observations. Though SCREAM cloud fraction in Figure 11a may appear to be underrepresented, we note that its deficiencies are small compared to most GCMs (Bogenschutz et al., 2021). In addition, cloud liquid water in Figure 11b matches observations almost perfectly. Figure 11c depicts a snapshot of the SWCRE on March 1, 2020 at 18:00:00 UTC from SCREAM to demonstrate the model's ability to simulate healthy coastal Sc cloud decks and the gradual transition to more broken cloud.
Enlarge this image.Figure 11. Temporally and spatially averaged profiles of (a) cloud fraction and (b) cloud liquid water for SCREAM and C3M. SCREAM profiles are averaged for the month of February 2020 while the C3M represents the February climatology from 2006 to 2010. Both SCREAM and C3M profiles represent spatial averages from the southeast Pacific coastal stratocumulus region bounded from 35° to 15°S and 275° to 290°E. The area used for spatial averaging is denoted in (c), which represents a snapshot of shortwave cloud radiative effect from SCREAM for March 1, 2020 at 18:00:00 UTC.
Figure 12 displays the temporally averaged curtain of cloud fraction along the 20°S transect across the stratocumulus-to-deep-convection transition for SCREAM February 2020 average and C3M February climatology from 2006 to 2010. When read from right to left (i.e., along the direction of prevailing easterly winds), C3M observations depict a gradual deepening of cloud in the lower troposphere over progressively warmer SSTs. SCREAM looks reasonable near the coast but fails to deepen to the W and is generally too thin in depth and too weak. This was an unintended consequence of tuning choices made in the SHOC parameterization to achieve reasonable radiation balance; further tuning since this simulation has improved the realism of trade cumulus.
Enlarge this image.Figure 12. Temporally averaged curtain of cloud fraction along the 20°S transect across the stratocumulus to deep convection transition. SCREAM clouds are averaged over the month of February 2020 while the C3M represents the February climatology from 2006 to 2010. Both SCREAM and C3M profiles represent curtains bounded from 24°S to 16°S.
Evaluating the spatial distribution of precipitation from a 40-day simulation is challenging. Forty days is too long for comparison against weather events but too short to average out the effects of individual storms. Zonal-averaging beats down some of this weather noise and large-scale tropical precipitation structure is probably robust, but results should still be taken with caution. In Figure 13, zonal-average precipitation is found to generally agree well with both GPM and ERA5 except for excessive rainfall on the equatorward side of the northern-hemisphere storm track and at the poleward edges of the tropics. GPM is known to be biased low at higher latitudes due to problems detecting light rain and snow (G. Huffman et al., 2019), which might partially explain storm track and polar biases. Figure 14 shows that tropical zonal-mean bias is due to a complicated mixture of differences in the meridional structure of precipitation. SCREAM tends to have stronger precipitation on the east side of land masses, in particular over the Maritime Continent (which has been a long-standing bias in E3SM; Golaz et al., 2019) and west of Madagascar. Heavy precipitation in the ITCZ extends too far east, which is another persistent E3SM bias. Precipitation in the South Pacific Convergence Zone (SPCZ) is, on the other hand, too weak and a bit too zonal. This may indicate that SCREAM (like most climate models) suffers from double-ITCZ problems (Li & Xie, 2014). Precipitation over the Amazon rain forest is slightly too strong, which is the opposite from what is seen in conventional climate models (Yin et al., 2012).
Enlarge this image.Figure 13. Zonal-average precipitation averaged over the last 30 days of the SCREAM simulation.
Enlarge this image.Figure 14. Tropical precipitation over the last 30 days of the SCREAM run (top), GPM observations averaged over the same period (middle), and their difference (SCREAM minus GPM, bottom).
A great success of cloud-resolving models are their ability to simulate the diurnal cycle of precipitation (Khairoutdinov et al., 2005; Sato et al., 2009; Stevens et al., 2019). This is a feature which coarser resolution models struggle with (Covey et al., 2016), though progress has been reported (Xie et al., 2019). As documented in Figure 15, SCREAM is able to capture the morning-time peak over the oceans and late afternoon peak over land. The diurnal cycle over the Maritime Continent and Madagascar—two areas dominated by sea breezes—is actually stronger than observed in GPM (but is weaker in magnitude than TRMM's observed climatology; not shown). Stronger diurnal amplitude in these areas is perhaps unsurprising given that daily mean precipitation was also noted to be too high in these regions.
Enlarge this image.Figure 15. Diurnal cycle of precipitation averaged over the last 30 days of the SCREAM run (top) and GPM observations (bottom). Hue indicates time of peak precipitation and intensity indicates diurnal amplitude. Amplitudes less than 1 mm day−1 are colored white and colors saturate at 25 mm day−1.
Like conventional GCMs (Na et al., 2020; Stephens et al., 2010), SCREAM has a tendency toward having too much drizzle and not enough strong precipitation (Figure 16). The magnitude of this bias is, however, much smaller than typically found in conventional GCMs (e.g., Caldwell et al., 2019; Figure 12). Thus, we consider simulation of heavy precipitation to be a victory for SCREAM.
Enlarge this image.Figure 16. Histogram of tropical precipitation over the ocean for the last 30 days of the SCREAM simulation compared to equivalent days from GPM.
Hovmöller diagrams showing precipitation averaged from 5°N to 5°S latitude as a function of longitude and time are useful for evaluating the temporal intermittency and propagation of tropical convection which collectively result in the MJO (Madden & Julian, 1971). Usually MJO analyses filter out signals outside of a 20–90 days window, but our 40-day simulation precludes such processing. A longer simulation is needed for statistical robustness, but it seems clear in Figure 17 that SCREAM triggers convection too frequently. This feature is also apparent in instantaneous snapshots of precipitation, water vapor, and cloud mass (not shown). We are still investigating the source of this “popcorn convection”, which also appears in other convection-permitting regional and global models (Arnold et al., 2020; Kendon et al., 2012). As found for other GCPMs (Miura et al., 2007; Miyakawa et al., 2014), SCREAM does a good job of propagating convective events eastward, though its propagation speed is perhaps slightly fast.
Enlarge this image.Figure 17. Precipitation averaged from 5°N to 5°S as a function of longitude (x-axis) and time (y axis) from SCREAM (left) and GPM precipitation observations (right).
The statistical analysis of precipitation above is important, but it ignores the fact that precipitation comes from storms whose characteristics vary regionally. The next few subsections explore SCREAM's treatment of important storm types.
Tropical CyclonesTropical cyclones (TCs) are some of the most intense storms in the world, combining intense precipitation with winds frequently in excess of 30 m s−1. Although some global models are able to represent TC frequency and intensity well at 0.25° grid spacing, resolving the inner structure of these storms requires much finer resolution (Judt et al., 2021; Wehner et al., 2014; Zarzycki & Jablonowski, 2015). A key advantage of running global convection-permitting models is the ability to represent and study multiscale interactions between the inner structure of tropical cyclones and the large-scale environment (Satoh et al., 2019). In the first phase of the DYAMOND project, models produced a wide range of tropical cyclone counts and intensities with counts as low as four to as high as 20, while in reality there were 14 (Judt et al., 2021; Stevens et al., 2019). In this section, we provide a brief and broad overview of the tropical cyclones identified in the SCREAM simulation.
Figure 18 shows TC tracks during the simulation period from SCREAM, ERA5, and IBTrACS observations (Knapp et al., 2010, 2018). SCREAM and ERA5 tracks are computed using the TempestExtremes (TE) algorithm and the criteria described in Appendix A1, while IBTrACS are based on expert judgment. Large discrepancies between IBTrACS and reanalysis data sets highlight the importance of using consistent criteria to classify storms. Note as well that the chaotic nature of weather means that storms later in the SCREAM simulation are not expected to match those found in ERA5 or in observations. Within the days of potential predictability (up to 2 weeks), one TC exists in both the SCREAM simulation and ERA5 data (Moderate Tropical Storm Diane). Another storm that is present in ERA5 (Moderate Tropical Storm Esami) does not organize in SCREAM, although a weak low pressure region does persist. Over the entire simulation period, we identify five tropical cyclone tracks in SCREAM during the 40-day simulation and six tracks in the ERA5 reanalysis data. All five TCs in SCREAM occur in the Southern Hemisphere, with four over the Indian Ocean and one off the northwestern coast of Australia over the Pacific Ocean, all broadly located where TCs are found in the reanalysis.
Enlarge this image.Figure 18. Tracks of tropical storms from IBTrACS (gray) and identified by the Tempest Extremes algorithm in SCREAM (purple) and in ERA5 (orange) between January 20 and February 28, 2020. Starting location is indicated with a plus (+).
In ERA5, Diane starts off as a tropical depression with central pressure of 1,020 hPa, but its pressure drops down to 990 hPa by January 26 with sustained maximum winds of 25 m/s (49 knots) or more. The simulated storm track from SCREAM closely follows that found in the reanalysis (Figure 19a), although it forms farther to the east and moves eastward more slowly. The maximum wind speed within a 6° × 6° box around the storm is also higher in the model, but this is likely due to the use of native grid data in SCREAM and the coarser regridding of the reanalysis data. Precipitation rates in SCREAM and reanalysis closely follow each other until ERA5 starts tapering off while SCREAM continues growing. GPM precipitation, however, includes a period of much stronger precipitation which isn't captured by either model simulation. Interestingly, SCREAM has a strong and regular diurnal cycle of precipitation which isn't found in the other timeseries. Although the data for this storm from IBTrACS spans a much shorter period in the storm lifetime than identified by Tempest Extremes from either the reanalysis or the SCREAM simulation, the magnitudes of central pressure and maximum 10-m wind speed are in good agreement between SCREAM, ERA5, and IBTrACS for the period that does overlap.
Enlarge this image.Figure 19. (a) Tracks of the tropical storm Diane from IBTrACS (gray) and as identified by the Tempest Extremes algorithm in SCREAM (purple) and in ERA5 (orange). Starting location is indicated with a + symbol. Shown below the tracks are time evolution of the storm's minimum central pressure (b), maximum 10-m wind speeds within 5° of the storm center (c), and area-averaged precipitation rate (d). (e–h) Same as (a–d) but for February 10 tropical cyclone in SCREAM simulation. No observational equivalent is shown, because it is outside the period of predictability.
Because Severe Tropical Storm Diane does not fully develop a canonical tropical cyclone structure and exhibits hurricane force winds only for a few hours, we take a more detailed look at a stronger storm in the model which forms on Feb 10 and produces surface wind speeds which classify it as a category 3 hurricane (Figure 19g). For reference, the storm's maximum intensity (based on minimum surface pressure values) is the median of the five storms tracked in SCREAM (not shown). Figure 19e shows the cyclone track, which spans 16 days. The surface pressure rapidly drops from February 11 to February 14, a minimum pressure of 930 hPa on Feb 16, when maximum 10-m wind speeds are also reached. By that point, the storm has formed a distinctive eye, ringed by strong precipitation rates reaching 100 mm/hr and wind speeds greater than 60 m/s (Figure 20). The high surface wind speeds drive surface latent heat fluxes greater than 500 W m−2, and a vertical north-south curtain centered on the point of minimum surface pressure shows the boundary layer flow is transporting energy toward the eye, particularly in the southern half of the storm (Figure 20).
Enlarge this image.Figure 20. Instantaneous planar and curtain view of February 11 tropical cyclone at maximum intensity on February 16 0 UTC. On the left column are planar views of the (a) outgoing longwave radiation, (b) precipitation rate, (c) latent heat flux, and (d) 10-m wind speed. On the right column is a north-south curtain snapshot through the center of the storm of the (e) cloud liquid mass mixing ratio, (f) ice mass mixing ratio, (g) zonal wind speed, and (h) meridional wind speed.
More analysis is necessary for an in depth study of the storm characteristics in SCREAM, as was done by Judt et al. (2021) for the models participating in the first phase of DYAMOND. However, as Figures 19 and 20 indicate, SCREAM produces tropical cyclones with reasonable eye-wall structure and adequate surface wind intensities, which provide promise for future attempts to simulate observed tropical cyclones using the model.
Extratropical CyclonesIn mid- and high-latitude regions, extratropical cyclones (ECs) are a large source of day-to-day weather variability. ECs are a major pathway for water evaporated from the ocean to precipitate over land; Hawcroft et al. (2012) suggest that as much as 90% of the surface precipitation along midlatitude storm tracks is attributed to ECs. ECs are also behind a majority of extreme precipitation events, particularly in the northeast US where ECs are responsible for more than 80% of winter-time extreme precipitation (Agel et al., 2015; Pfahl & Wernli, 2012). With increasing resolution, ECs are better represented in global models (Jung et al., 2006), and a recent study using a set of global storm-resolving model simulations shows an increase of 7%/K in precipitation rate from the most intense extratropical cyclones with warming, which differs from the 2%–3%/K increase expected in the global mean (Kodama et al., 2019).
Over the simulation time period, 87 ECs are identified in SCREAM and 80 are found in ERA5 using the TempestExtremes algorithm (see Appendix A2 for details). Their geographic distributions in the Southern and Northern Hemispheres are shown in Figures 21a and 21b. In the Northern Hemisphere, the density of storms in both SCREAM and ERA5 is largest over the Atlantic and Pacific Ocean basins, with many storms originating close to the western boundary currents. This is consistent with observed climatologies of cyclone statistics (Sinclair, 1997). Bomb cyclones (ECs with surface low pressures dropping more than 24 hPa over a 24 hr period [Sanders & Gyakum, 1980]) are present in both SCREAM (11) and in ERA5 (15). While small numbers prevent us from making conclusive statements, spatial distributions in ERA5 and SCREAM seem consistent.
Enlarge this image.Figure 21. Geographic distribution of extratropical cyclones identified in SCREAM (cyan) and ERA5 (orange) using the TE algorithm (described in Appendix A2) for the Southern Hemisphere (a) and Northern Hemisphere (b). Dark blue tracks indicate bomb cyclones in SCREAM, whereas brown tracks indicate bomb cyclones in ERA5. (c) The latitudinal distribution of 6 hourly snapshots of extratropical cyclones in ERA5 (black) and SCREAM (red). The dashed black line indicates the distribution found in ERA5 for the DYAMOND2 period (January 20 through February 28, 2020). Solid black line indicates the average distribution for January 20 to February 28 of 1979 through 2020 in ERA5 with gray shading indicating maximum and minimum ranges for each year.
Figure 21c shows the frequency of ECs by latitude band. ECs are counted separately in each 6 hourly snapshot in this plot, so counts in this plot are much higher than the ∼80 storms quoted above for SCREAM and ERA5, which tracked single storms across time. In both hemispheres, SCREAM has a more peaked distribution with maximum frequency at the upper limit of the observed count from the 1979–2020 period. The excessively peaked EC count structure in the northern hemisphere is consistent with zonal precipitation bias shown in Figure 13. Interestingly, modeled southern hemisphere storm track precipitation in Figure 13 matches ERA5 almost perfectly despite having excessive EC count around 50°S. Storm composites show that Southern Hemisphere extratropical cyclones in SCREAM are associated with less rain than ERA5, which might explain this apparent paradox (not shown). Peak latitude is roughly consistent with observations in each hemisphere, though is displaced slightly poleward in the northern hemisphere.
We noted earlier that large swaths of the Southern Ocean in SCREAM have too much absorbed shortwave radiation compared to CERES-EBAF retrievals (Figure 9). Many climate models share biases where the cold sector of storms does not reflect enough incoming shortwave radiation, while the warm sector is less biased (Bodas-Salcedo et al., 2014). To examine whether this is the case in SCREAM, we construct composites of the cyclones tracked in SCREAM between 40° and 60°S. This latitude band is consistent with those of Bodas-Salcedo et al. (2014), but ignores storms with centers poleward of 60°S (to remove complications due to the reflectivity of sea ice). Figure 22 shows the composite of the pseudo-cloud albedo for SCREAM and its difference with CERES-SYN-based estimates. The pseudo-cloud albedo is defined here as the shortwave cloud radiative effect divided by the local solar insolation. By using a pseudo-cloud albedo rather than reflected shortwave radiation, we remove the potential impact of biases in the latitudinal distribution of ECs on our assessment of SCREAM's cloud reflectivity. Indeed, like the GCMs studied by Bodas-Salcedo et al. (2014), there is less cloud reflection in the cold sector of SCREAM's storms (−4.9% in the cold western half of the storm), compared to the storms captured in ERA5. However, the warm-sector of the storm also shows lower cloud albedo (−3.8% in the warm eastern half of the storm), showing that in SCREAM, there is a general lack of cloud reflection. Fixing this bias is a research priority.
Enlarge this image.Figure 22. (left) Composite of the pseudo-cloud albedo in extratropical cyclones found between 45° and 60°S in SCREAM (January 20–February 28). Composites average over all 6 hourly snapshots centered on identified ECs and are plotted such that north is oriented upward. (right) Difference in storm composite pseudo-cloud albedo between storms in SCREAM and in reanalysis and satellite data (ERA5/CERES-SYN) for storms occurring during the same period.
Atmospheric rivers (ARs) are long, narrow, and transient corridors of enhanced vapor transport typically associated with the low-level jet stream ahead of the cold front of an extratropical cyclone (AMS, 2019). As noted by Zhu and Newell (1998), ARs are responsible for approximately 90% of poleward vapor transport. Water resources in the western U.S. are strongly tied to ARs, with landfalling ARs providing approximately 20%–50% of total wet season precipitation (Dettinger et al., 2011; Lavers & Villarini, 2015) and 30%–40% of mountain snowpack (Guan et al., 2010). One such landfalling atmospheric river observed in the SCREAM simulation along the west coast of North America is depicted in Figure 23.
Enlarge this image.Figure 23. Snapshot of a landfalling atmospheric river along the west coast of North America that occurs on February 11, 23:00:00 UTC. Vertically integrated water vapor is indicated in transparent grayscale with opaque/white regions having integrated vapor greater than 40 kg m−2. Colors indicate precipitation intensity.
To assess the quality of ARs in the SCREAM simulation, we track ARs over the simulation period using the TempestExtremes atmospheric river detection and tracking algorithm (McClenny et al., 2020; Ullrich & Zarzycki, 2017) as described in Appendix A3. In Figure 24 the properties of these tracked features are then compared to analogously tracked features from all January 20 through Mar 28 periods in ERA5 data (1979–2020), roughly following the approach discussed in Rutz et al. (2019). In general SCREAM falls well within the climatological range from ERA5 historical simulations, except for a slight underestimation of AR frequency around 35° north and south of 50°S. For 2020, ERA5 predicts abnormally high AR activity while SCREAM is slightly weaker than ERA5's long-term average. Without an ensemble of simulations to compare against, however, such a discrepancy could very easily be attributed to interannual variability.
Enlarge this image.Figure 24. Properties of tracked atmospheric rivers in both hemispheres over the period January 20 through February 28 of each year in (red) the SCREAM DYAMOND2 simulation and (gray shaded region with mean shown with black solid line) 1979–2020 ERA5 reanalysis. Plots refer to (top) average atmospheric river frequency, as a percent of the full longitudinal band, with results from 2020 depicted with a black dashed line; (middle) zonally averaged northward integrated vapor transport (IVTn) at grid points flagged as part of/not part of atmospheric rivers; (bottom) mean fractional contribution of northward vapor transport from atmospheric rivers relative to all northward vapor transport.
The underestimation of AR frequency in southern high latitudes is associated with anomalously low eastward integrated vapor transport (IVT), which is in turn due to anomalously low eastward wind speeds compared to ERA5 (not shown). Interestingly, Figure 21 shows that EC frequency was actually too high where we find AR frequency to be too low. Perhaps ECs are spending too much time in this region due to low wind speeds? Nonetheless, the fractional contribution of ARs to poleward transport of moisture is almost identical to the climatological mean performance from ERA5, suggesting consistency of the underlying physical processes. Overall we conclude that SCREAM performs well in its representation of ARs and their associated contribution to poleward transport of vapor.
Cold-Air OutbreaksMarine cold air outbreaks (MCAOs) occur when cold air of polar or continental origin flow over warm ocean waters. Because of the strong air-sea temperature differences and typical higher surface wind speeds, cold air outbreaks are regions of strong surface turbulent heat fluxes that can reach 1,000 W m−2 (Shapiro et al., 1987) and can impact frontogenesis (Terpstra et al., 2016). General circulation models (GCMs) have, however, not represented clouds under these conditions very well (Rémillard & Tselioudis, 2015). The models tend to simulate too little stratiform cloud cover in these regions (Bodas-Salcedo et al., 2014; Field et al., 2014). In this section, we describe the frequency and intensity of MCAOs in the SCREAM simulation relative to reanalysis (ERA5) during the same time period and examine the surface flux and cloud properties for a single cold air outbreak event that occurs early in the simulation over the Kuroshio current.
To identify and quantify cold air outbreaks, we use the cold air outbreak index (M) as described by Fletcher et al. (2016), which is quantified as the potential temperature difference between the surface skin and 800 hPa. Any oceanic region with a positive value of M denotes a region undergoing a cold air outbreak. If we compare the frequency of cold air outbreaks in SCREAM and in ERA5 over the global oceans, we see general agreement of where and how often cold air outbreaks occur (Figures 25a and 25c). Cold air outbreaks tend to occur most prominently in the winter Northern Hemisphere along the eastern edges of continents and southern edges of the sea-ice. In regions where SCREAM produces cold air outbreaks (e.g., over the Kuroshio current, Gulf stream current, and south of Alaska), M frequency tends to be higher. MCAOs are, however, greatly underestimated to the south and east of Greenland. This is unsurprising since 2-m temperature is far too warm over Greenland (Figure 6), likely due to meridional wind biases discussed in Section 6.2. Except for a slight overestimation, SCREAM also tends to capture well the intensity of the strongest of cold air outbreaks (Figures 25b and 25d).
![View Image - Figure 25. Frequency of cold air outbreaks (based on the M of Fletcher et al. [2016]) in SCREAM over the month of February 2020 in SCREAM (a) and in ERA5 (c). Also shown is the 95th percentile value of M (including non-cold air outbreak instances) over the same period in SCREAM (b) and ERA5 (d).](https://media.proquest.com/media/hms/THUM/34/hGXgL?_a=ChgyMDI2MDMyMTE5NDEyODMzMjoyMzU3NzESBzEzODMzMjMaCk9ORV9TRUFSQ0giDjIxNi43My4yMTYuMTM2KgY2MTY2NjcyCjI2MDA4MzM2NDU6C0lubGluZUltYWdlQgEwUgZPbmxpbmVaA0lMSWIEVEhVTWoKMjAyMS8xMS8wMXIKMjAyMS8xMS8zMHoAggEgUC0xMDAwMDg0LTIwMTM5NS1TVEFGRi1udWxsLW51bGySAQZPbmxpbmXKAWdNb3ppbGxhLzUuMCBBcHBsZVdlYktpdC81MzcuMzYgKEtIVE1MLCBsaWtlIEdlY2tvOyBjb21wYXRpYmxlOyBDbGF1ZGVCb3QvMS4wOyArY2xhdWRlYm90QGFudGhyb3BpYy5jb20p0gESU2Nob2xhcmx5IEpvdXJuYWxzmgIHUHJlUGFpZKoCG09TOkVNUy1GdWxsVGV4dC1nZXRGdWxsVGV4dMoCD0FydGljbGV8RmVhdHVyZeICAPICAPoCAVmCAwNXZWKKAxxDSUQ6MjAyNjAzMjExOTQxMjgzMzI6NjE3NjY1&_s=YvpzxcGD3bk53nxazVWcfrRGuYU%3D "View Image - Figure 25. Frequency of cold air outbreaks (based on the M of Fletcher et al. [2016]) in SCREAM over the month of February 2020 in SCREAM (a) and in ERA5 (c). Also shown is the 95th percentile value of M (including non-cold air outbreak instances) over the same period in SCREAM (b) and ERA5 (d).")
Enlarge this image.Figure 25. Frequency of cold air outbreaks (based on the M of Fletcher et al. [2016]) in SCREAM over the month of February 2020 in SCREAM (a) and in ERA5 (c). Also shown is the 95th percentile value of M (including non-cold air outbreak instances) over the same period in SCREAM (b) and ERA5 (d).
To study the cloud fields that form under the simulated cold air outbreaks in SCREAM, we focus on a cold air outbreak event that flows off the Asian continent over the Kuroshio current from January 21 to January 22. We examine the cold air outbreak characteristics over the 24 hr period of January 22 to exclude any impacts of the cold front. The simulated sensible heat flux generally matches ERA5, but is a bit too smooth and too big (Figures 26a and 26d). Good spatial agreement may be an artifact of prescribed SST; smooth features are probably due to use of a coarser (∼6 km) ocean grid in this region. Excessive magnitude is unsurprising given surface wind speed biases mentioned in Section 6.2 and again apparent from comparing Figure 26 panels (b vs. e). Surface air temperature bias does not contribute to excessive surface fluxes (not shown).
Enlarge this image.Figure 26. The daily mean sensible heat flux over the Kuroshio region bounded from 29° to 49°N and 141.5° to 171.5°E in SCREAM (a) and ERA5 (d) for the cold air outbreak on January 22. Also shown are similar daily mean values of 10-m wind speed (b) SCREAM; (e) ERA5 and shortwave cloud radiative effect (c) SCREAM; (f) CERES-SYN.
Although GCMs tend to underestimate the occurrence of MCAO clouds and SCREAM itself was shown earlier to suffer from a deficiency in clouds in other regimes, a comparison of the shortwave cloud radiative effect between the model and CERES-SYN suggests good agreement in the MCAO regime (Figures 26c and 26f). In Figure 27 we take a closer look at cloud structure by comparing a snapshot of shortwave cloud radiative effect from SCREAM against a visible satellite image taken at the same time from Himawari-8 (Bessho et al., 2016). Similarity between the observed and simulated cloud structures is striking, particularly since this image is taken 2 days into a free-running simulation. In particular, cloud streets in SCREAM form along the direction of the flow, before transitioning into broken and open-cellular convection further offshore. The model's ability to capture this transition suggests that SCREAM's combination of resolution and boundary layer/cloud parameterizations contains the physics necessary to capture cloud transitions in cold air outbreaks. Further analyses compositing many cold air outbreak events would be necessary to draw more general conclusions.
Enlarge this image.Figure 27. Cold-air outbreak off Siberia on January 22, 2020 at 2:00:00 UTC (∼local noon) from a Himawari visible satellite image (left) and shortwave cloud radiative effect from SCREAM (right). Visualization is over a region bounded by 29°–49°N and 141.5° to 171.5°E.
The overall takeaway from this work is that 3.25 km global models solve a lot of the long-standing problems in global climate modeling even without the detailed optimization and tuning which is typically so important for GCM skill. In particular, SCREAM does an excellent job simulating precipitation; its diurnal cycle (Figure 15) and intensity distribution (Figure 16) are particularly realistic. Tropical and extratropical storm frequency and structure (Sections 6.5–6.7) are also impressive. The vertical structure of tropical convection (Figure 10) is also much improved relative to typical GCMs. Coastal stratocumulus (Figure 11) and cold-air outbreaks (25–27), which are perennially difficult to simulate not just in GCMs (Rémillard & Tselioudis, 2015) but also in limited-area CPMs (Klein et al., 2009), are also well captured. We suspect that the SHOC cloud/turbulence parameterization and fine vertical resolution within SCREAM were important for this success.
Several biases in SCREAM are familiar from conventional GCMs. Clarifying whether these biases are caused by processes unresolved at 3.25 km grid spacing would be a large step toward understanding and therefore fixing these perennial problems. One such bias is the tendency for the South Pacific Convergence Zone to be too zonal (Figures 13 and 14). This suggests that resolution doesn't resolve the double-ITCZ bias that plagues lower-resolution models. This finding is consistent with the result of Stevens et al. (2019) for other GCPMs. Another bias in lower-resolution versions of E3SM which persists in SCREAM is a tendency for precipitation in the West Pacific to be maximized over the Maritime Continent rather than to the east over the ocean.
Analysis for this paper also revealed several deficiencies which will be fixed in future model versions. First, cloud fraction near the tropopause is corrupted by the use of a relative-humidity based ice cloud fraction scheme tuned for low resolutions (Figure 4). Because these spurious clouds had no mass, they had little practical impact on the simulation, but users of SCREAM DYAMOND2 data should be careful to use our post-facto-generated cloud-mask-based cloud fraction for future analysis. Overly strong surface wind speed is a second deficiency (Figure 7). Upper level winds are generally reasonable but have unrealistic poleward transport south of Greenland and around Australia (Figure 8). Surface temperature at high latitudes is also problematic (Figure 6). One potential reason for this is a land initial condition with low snowpack in mountainous regions exacerbated by potentially poor tuning of the lower limit on turbulent mixing in stable conditions and aforementioned biases related to heat transport into polar regions. Another issue is a prevalence of frequent, small ″popcorn″ convective events (Figure 17). Finally, cloud tuning should be improved. Shortwave reflection and longwave emission are too weak (Figure 9) and low-level clouds tend too much toward stratus and too little toward shallow convection (Figure 12). Issues like these are expected for a new model version and many of these issues have an obvious solution. We are releasing this initial model without fixing these problems to match the timing of the DYAMOND2 intercomparison, because there will always be something else to fix, and because using a model for science and writing papers is by far the fastest way to find problems.
This simulation is a milestone rather than an endpoint in SCREAM development. In addition to fixing the issues identified above, the major focus of the SCREAM project is on completing the computationally performant C++ implementation of the model. We hope to perform longer, more realistic simulations soon.
For feature tracking in the DYAMOND2 simulation we use TempestExtremes 2.1 (Ullrich & Zarzycki, 2017), available from ZENODO at
Tropical cyclone tracking is performed on 6-hourly data following (Zarzycki et al., 2017). The search is performed for local minima in the sea level pressure (PSL) which are accompanied by an increase of 200 Pa over a distance of 5.5° great circle distance (GCD). Tropical cyclones are further defined by the presence of an upper-level warm core which is characterized by anomalous thickness in the geopotential height between 500 and 200 hPa. Here we require that this thickness drop by 6.0 m over a distance of 6.5° GCD, where the maxima in the layer thickness must be within 1.0° GCD of the pressure minima. Following this only the most intense features within 6.0° GCD are retained. Tracks are then stitched together in time, where sequential features must be within 8.0° GCD, must persist for at least 10 time steps (2.5 days), can have no more than 3 sequential 6-hourly time steps where no detection is found, must have a 10 m wind speed greater than 10 m s−1 for at least 10 steps along the trajectory, and must be within 50°S and 50°N for at least 10 steps along the trajectory. The commands are as follows:
-
$TEMPESTEXTREMESDIR/DetectNodes −−in_data_list DYAMOND_TC_files.txt
-
−−out DYAMOND_DN.txt −−searchbyminPSL
-
−−closedcontourcmd ″PSL, 200.0, 5.5, 0;_DIFF(Z200, Z500), −6.0, 6.5, 1.0″
-
−−mergedist6.0 −−outputcmd ″PSL, min, 0; WINDSPD_10M, max, 2″ −−timefilter ″6 hr″
-
$TEMPESTEXTREMESDIR/StitchNodes −−in DYAMOND_DN.txt
-
−−out DYAMOND_TC_tracks.txt −−in_fmt "lon,lat,slp,wind" −−range 8.0
-
−−mintime "10" −−maxgap "3"
-
−−threshold "wind, ≥, 10.0, 10; lat, ≤, 50.0, 10; lat, ≥, −50.0, 10"
As with tropical cyclones, extratropical cyclone tracking is performed on 6-hourly data. Candidates are first detected as minima in the difference between the sea-level pressure (PSL) and the average sea-level pressure over the entire simulation (PSL_climo). We require that this difference increase by 200 Pa within 5.5° GCD of the candidate. We further eliminate points that have an upper-level warm core, as these are likely tropical cyclones, by removing candidates with a drop in the 500–200 hPa layer thickness of 6.0 meters within 6.5° GCD of the point of maximum layer thickness within 1.0 degrees of the candidate. Following this only the most intense features within 6.0° GCD are retained. Tracks are then stitched together in time, where sequential features must be within 8.0° GCD, must persist for at least 8 time steps (2.0 days), can have no more than 2 sequential 6-hourly time steps where no detection is found, must have a surface geopotential less than 700.0 for at least 8 time steps, and must have a distance of 6.0° GCD between genesis and termination point. The commands for these operations are as follows:
-
$TEMPESTEXTREMESDIR/bin/DetectNodes −−in_data_listDYAMOND_ETC_files.txt
-
−−out DYAMOND_DN_ETCs.txt −−searchbymin ″_DIFF(PSL, PSL_climo)″ −−timefilter ″6 hr″
-
−−closedcontourcmd ″_DIFF(PSL, PSL_climo), 200.0, 5.5, 0″
-
−−noclosedcontourcmd ″_DIFF(Z200, Z500), −6.0, 6.5, 1.0″ −−mergedist 6.0
-
−−outputcmd ″PSL, min, 0; _DIFF(PSL, PSL_climo), min, 0; WINDSPD_10M, max, 5; PHIS, min, 0″
-
$TEMPESTEXTREMESDIR/bin/StitchNodes −−in DYAMOND_DN_ETCs.txt
-
−−out DYAMOND_ETC_tracks.txt −−in_fmt "lon,lat,psl,pslanom, wind,phis"
-
−−range 8.0
-
−−mintime "8" −−maxgap "2" −−min_endpoint_dist 6.0 −−threshold "phis, ≤, 700, 8"
Atmospheric river tracking is performed using the tracker employed in McClenny et al. (2020). Grid points poleward of 15° N/S are flagged where the Laplacian of the integrated vapor transport (evaluated using eight points with radius 10° GCD) is less than 20,000 kg m−1 s−1 rad−2. Only contiguous regions with area greater than 4 × 105 km2 are retained in this operation. Since high IVT blobs can include tropical cyclones, we also remove all points within 10° GCD of TCs detected using the method described in Appendix A1. The commands for these operations are as follows:
-
$TEMPESTEXTREMESDIR/DetectBlobs −−in_data CAT_TUQ, TVQ_256×512.eam.nc
-
−−outCAT_ARs_256×512.eam.nc−−minabslat 15 −−geofiltercmd ″area, >=, 4e5km2″
-
−−thresholdcmd ″_LAPLACIAN{8,10} (_VECMAG(TUQ, TVQ)),<=,−20,000, 0″
-
$TEMPESTEXTREMESDIR/NodeFileFilter −−in_nodefile DYAMOND_TC_tracks.txt
-
−−in_fmt "lon,lat" −−in_data CAT_ARtag_256 × 512.eam.nc
-
−−out_data CAT_ARtag_TCfiltered_256×512.eam.nc −−var "binary_tag"
-
−−bydist 10.0 −−invert
This research was supported as part of the Energy Exascale Earth System Model (E3SM) project, funded by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research. It used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231. This research also used resources from the Argonne Leadership Computing Facility at Argonne National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under contract DE-AC02-06CH11357. This paper was prepared by LLNL under Contract DE-AC52-07NA27344. H.-Y. Ma is funded by the Regional and Global Model Analysis program area and Atmospheric System Research program of the U.S. Department of Energy. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. We acknowledge and thank Jingyu Wang for accessing, formatting, and making the data available to us. The authors would also like to thank Jiwen Fan and Kai Zhang for their contributions to fixing bugs in our P3 microphysics implementation and three anonymous reviewers for excellent suggestions to improve the manuscript.
Data Availability StatementThe SCREAM output described in this paper is publicly available as part of the DYAMOND2 intercomparison as described at
© 2021. This work is published under http://creativecommons.org/licenses/by/4.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.