Advancing precipitation prediction using a new-generation storm-resolving model framework – SIMA-MPAS (V1.0): a case study over the western United States

Global climate models (GCMs) have advanced in many ways as computing power has allowed more complexity and finer resolutions. As GCMs reach storm-resolving scales, they need to be able to produce realistic precipitation intensity, duration, and frequency at fine scales with consideration of scale-aware parameterization. This study uses a state-of-the-art storm-resolving GCM with a nonhydrostatic dynamical core – the Model for Prediction Across Scales (MPAS), incorporated in the atmospheric component (Community Atmosphere Model, CAM) of the open-source Community Earth System Model (CESM), within the System for Integrated Modeling of the Atmosphere (SIMA) framework (referred to as SIMA-MPAS). At uniform coarse (here, at 120 km) grid resolution, the SIMA-MPAS configuration is comparable to the standard hydrostatic CESM (with a finite-volume (FV) dynamical core) with reasonable energy and mass conservation on climatological timescales. With the comparable energy and mass balance performance between CAM-FV (workhorse dynamical core) and SIMA-MPAS (newly developed dynamical core), it gives confidence in SIMA-MPAS's applications at a finer resolution. To evaluate this, we focus on how the SIMA-MPAS model performs when reaching a storm-resolving scale at 3 km. To do this efficiently, we compose a case study using a SIMA-MPAS variable-resolution configuration with a refined mesh of 3 km covering the western USA and 60 km over the rest of the globe. We evaluated the model performance using satellite and station-based gridded observations with comparison to a traditional regional climate model (WRF, the Weather Research and Forecasting model). Our results show realistic representations of precipitation over the refined complex terrains temporally and spatially. Along with much improved near-surface temperature, realistic topography, and land–air interactions, we also demonstrate significantly enhanced snowpack distributions. This work illustrates that the global SIMA-MPAS at storm-resolving resolution can produce much more realistic regional climate variability, fine-scale features, and extremes to advance both climate and weather studies. This next-generation storm-resolving model could ultimately bridge large-scale forcing constraints and better inform climate impacts and weather predictions across scales.

1 Introduction

Climate models have advanced in many ways in the last decade including their atmospheric dynamical core and parameterization components. Advances in computer power have now enabled climate models to be run with nonhydrostatic dynamical cores at “storm-resolving” scales, on the order of a few kilometers (Satoh et al., 2019). These GSRMs (global storm-resolving models) have been constructed at a number of modeling centers (Satoh et al., 2019; Stevens et al., 2019; Dueben et al., 2020; Stevens et al., 2020; Caldwell et al., 2021). We expect an emerging trend in improving and applying the new modeling structures for a better and more accurate understanding of global and regional climate studies and weather-scale predictions.

The Community Earth System Model (CESM) has been used in a wide range of climate studies. For high-resolution CESM applications (but hydrostatic only), the variable-resolution (VR) CESM-SE (spectral element core) for regional climate modeling has been used in many regional climate studies (such as Small et al., 2014; Zarzycki and Jablonowski, 2014, 2015; Rhoades et al., 2016; Huang et al., 2016; Huang and Ullrich, 2017; Bacmeister et al., 2018; Gettelman et al., 2018, 2019; van Kampenhout et al., 2019). Specifically, Rhoades et al. (2016) found that the VR CESM framework (with refinement at 0.25 and 0.125^∘ resolutions) can provide much enhanced representation of snowpack properties relative to widely used GCMs (such as CESM-FV 1^∘ and CESM-FV 0.25^∘) over the California region. Gettelman et al. (2018) found that the variable-resolution CESM-SE simulation (at 0.25^∘, ∼25 km) can produce precipitation intensities similar to the high-resolution simulation and has higher extreme precipitation frequency than the low-resolution simulation over the continental United States (CONUS) refinement region, close to observations.

More recently for storm-resolving model development, there have been two efforts to bring the dynamical core from the Model for Prediction Across Scales (MPAS) into CESM. The first effort involved implementing the hydrostatic atmospheric dynamical core in MPAS Version 1 in the Community Atmosphere Model (CAM), which is the atmospheric component of CESM. This effort made available the horizontal variable-resolution mesh capability of the MPAS spherical centroidal Voronoi mesh (Ringler et al., 2010) and led to a number of studies (e.g., Rauscher et al., 2013; Rauscher and Ringler, 2014; Sakaguchi et al., 2016). For example, Rauscher et al. (2013) found that tropical precipitation increases with increasing resolution in the CAM-MPAS using aquaplanet simulations.

Later, the static port of MPAS to CAM was updated with the nonhydrostatic MPAS atmospheric solver (Skamarock et al., 2012, 2014) to provide nonhydrostatic GSRM capabilities to CAM (Zhao et al., 2016). Neither of these ports was formally released, and the nonhydrostatic MPAS was not energetically consistent with CAM physics or its energy fixer given, among other things, the height vertical coordinate used by MPAS. Furthermore, the MPAS modeling system and its dynamical core, being separate from CESM, have evolved from these earlier ports. To address the issues in the earlier MPAS dynamical core ports to CAM/CESM, the MPAS nonhydrostatic dynamical core has been brought into CAM/CESM as an external component – i.e., it is pulled from the MPAS development repository when CAM is built – and all advances in MPAS are immediately available to CESM-based configurations using MPAS. This latest port was accomplished as part of the SIMA (System for Integrated Modeling of the Atmosphere) project. Importantly, this implementation also includes an energetically consistent configuration of MPAS, with its height vertical coordinate, the CAM hydrostatic-pressure coordinate physics, and the CAM energy fixer.

The MPAS dynamical core solves the fully compressible nonhydrostatic equations of motion and continues to be developed and used in many studies (Feng et al., 2021; Lin et al., 2022; also see https://mpas-dev.github.io/atmosphere/atmosphere.html, last access: April 2022). In this work, we test the storm-resolving capabilities in this new atmospheric simulation system. We use SIMA capabilities to configure a version of CESM with the MPAS nonhydrostatic dynamical core, called SIMA-MPAS instead of CESM-MPAS, since it is coupled only to a land model, with the other climate-system components being data components. In particular, we would like to answer the following question: can a nonhydrostatic dynamical core (dycore) coupled global climate model reproduce observed wet-season precipitation over targeted refinement regions? In addition, will this new development and modeling framework perform better or worse than a mesoscale model at similar resolution?

We aim to understand how this new SIMA-MPAS model configuration performs when configured for the storm-resolving (convection-permitting) scale for precipitation prediction over the western United States (WUS). Leveraging the recent significant progress in SIMA-MPAS development, we have undertaken experiments to understand the performance of SIMA-MPAS in precipitation simulations involving heavy storm events and relevant hydroclimate features at fine scales. We also explore large-scale dynamics and moisture flux transport over the subtropical region across the North Pacific. We evaluate the model results compared to both observations and a regional climate model. Employing the recent modeling developments in CESM with the MPAS dycore, the ultimate goal of this study is to evaluate the potential improvements to our understanding of atmospheric processes and predictions made possible with GSRM capabilities. We begin in Sect. 2 with a description of the model configurations and experiments. Section 3 describes the main results, including mean climatology diagnostics, precipitation and snowpack statistics, and large-scale moisture flux and dynamics. A summary and discussion follow in Sect. 4.

2 Methods, experiments, and dataset

2.1 Methods and experiments

As briefly mentioned in the Introduction section, we configure CESM2 (Danabasoglu et al., 2020) with the MPAS nonhydrostatic dynamical core and CAM6 physics. We call this configuration SIMA-MPAS. SIMA is a flexible system for configuring atmospheric models inside of an Earth system model for climate, weather, chemistry, and geospace applications (https://sima.ucar.edu, last access: April 2022). The components of this particular configuration also include the Community Land Model Version 5.0 (CLM5; with the MOSART river model) and prescribed observation-based SST (sea surface temperature) and ice. MPAS-Atmosphere employs a horizontal unstructured centroidal Voronoi tessellation (CVT) with C-grid staggering (Ringler et al., 2010), and its numerics exactly conserve mass and scalar mass. Both horizontal uniform meshes and variable-resolution meshes with smooth resolution transitions are available for MPAS-Atmosphere, and this study employs both mesh types. It uses a hybrid terrain-following height coordinate (Klemp, 2011).

We summarize here the key developments to the coupling of the MPAS dynamical core to CAM physics and changes to CAM physics to accommodate MPAS. Most of all, we would like to point out that a consistent coupling of the MPAS dynamical core with the CAM physics package is not trivial for several reasons. (1) MPAS uses a height (z)-based vertical coordinate, whereas CAM physics uses pressure. (2) The CAM physics package enforces energy conservation by requiring each parameterization to have a closed energy budget under the constant-pressure assumption (Lauritzen et al., 2022). For the physics–dynamics coupling to be energy consistent (i.e., not be a spurious source/sink of energy) requires the energy increments in physics to match the energy increments in the dynamical core when adding the physics tendencies to the dynamics state. When “mixing” two vertical coordinates, that becomes non-trivial. (3) The prognostic state in MPAS is based on a modified potential temperature, density, winds, and dry mixing ratios, whereas CAM uses temperature, pressure, winds, and moist mixing ratios for the water species. The conversion between (discrete) prognostic states should not be a spurious source/sink of energy either. (4) Lastly, the energy fixer in CAM that restores energy conservation due to updating pressure (based on water leaving/entering the column), as well as energy dissipation in the dynamical core and physics–dynamics coupling errors (see Lauritzen and Williamson, 2019), assumes a constant-pressure upper boundary condition. MPAS assumes a constant height at the model top, so the energy fixer needs to use an energy formula consistent with the constant volume assumption. The details of the energy-consistent physics–dynamics coupling and extensive modifications to CAM physics to accommodate MPAS are beyond the scope of this paper and will be documented in a separate source.

In terms of scale awareness, there are two aspects related to the model physics in the configuration that must be considered when employing regionally refined meshes. First, features resolvable in the finer regions of the mesh may not be resolvable in the coarser regions of the mesh. These features, e.g., deep convection in this study, need to be parameterized in the coarse mesh regions and not parameterized in the fine mesh regions, typically with the parameterization reducing its adjustment gradually in the mesh transition regions. Second, the time step used for the physics is the same over the entire mesh. i.e., in both coarse and fine regions, and the time step in CESM-MPAS is chosen to be appropriate for the smallest grid, as indicated in Table 1. Within our simulations, the balance of deep convective (diagnostic) and stratiform (large-scale) precipitation changes with the mesh spacing. In addition, since the deep convective parameterization in CESM-MPAS has a closure with a fixed timescale, the parameterized convection produces less condensation in the coarse mesh regions compared to simulations with a larger time step appropriate for the coarser mesh (Gettelman et al., 2019). But in the simulations herein, most of the precipitation is strongly forced by the large-scale flow, with the larger condensation hypothesized to lead to larger rain rates. This is particularly important over the WUS complex terrains. The large-scale condensation scheme, part of the unified turbulence scheme (Golaz et al., 2002), has internal length scales that should adjust its distributions as the scale changes (less variance in the probability density functions, PDFs). Land-surface-related feedback is also resolution dependent with scale-aware surface heterogeneity and coupled land–atmosphere interactions to affect the phase and hydrological impacts resulting from the regional precipitation statistics.

With the above significant progress in SIMA-MPAS development, we would like to diagnose the performance of this new-generation model when applied at convection-permitting resolutions and when bridging both weather- and climate-scale simulations in a single global model. We have chosen the WUS (due to its hydroclimate vulnerability and complexity, heavily impacted by precipitation variability) as our study region to examine the precipitation features in SIMA-MPAS at fine scales during wet seasons. We aim to figure out when the model outperforms and underperforms when compared with a traditional regional climate model against the best available observations and observationally based gridded products at similar resolutions for mean and extreme precipitation. As mentioned in the Introduction, we would like to figure out whether a nonhydrostatic dycore coupled global climate model can reproduce observed wet-season precipitation over targeted refinement regions with heavy impacts. And will this new development and modeling framework perform better or worse than a mesoscale model at similar resolution? Those are important questions to answer given the long-standing biases in traditional hydrostatic GCMs for simulating heavy precipitation and extremes.

To answer those questions, we have designed and conducted a set of experiments as shown in Table 1. In detail, these can be described as follows.

• Set A. We have tested CESM2 at the same coarse resolution using both MPAS (at 120 km) as the nonhydrostatic core and finite volume (Danabasoglu et al., 2020) (at ∼1^∘) as the hydrostatic core for multiple years of climatology to obtain 5-year mean F2000 climatology (F2000climo, in which the SST and ice condition are prescribed at the same yearly climatology with a mean from the time period 1995–2005) at ∼1^∘ for both MPAS and the FV (finite-volume) dycore.

• Set B. As the main focus for this work, a variable-resolution mesh is configured with 3 km refinement centered over the WUS as shown in Fig. 1, for five wet-season simulations with 60–3 km mesh (years 1999 to 2004; mid-November to mid-March; FHIST component set for historical forcings, https://www.cesm.ucar.edu/models/cesm2/config/compsets.html, last access: April 2022); atmosphere conditions are initialized by Climate Forecast System Reanalysis (CFSR) data.

• Set C. In addition, we have also configured uniform 60 km simulations for two wet seasons in contrast to the 60–3 km ones (years 2000 to 2002, November to March).

• Set D. Lastly, to accommodate the recent changes to the Morrison and Gettelman (MG) microphysics scheme, we have also conducted simulations at 60–3 km resolution for the three wet seasons (years 1999–2002) using MG3 with graupel (Gettelman et al., 2019) instead of MG2 (Gettelman et al., 2015) as in the Set B simulations. Specifically, Gettelman et al. (2019; i.e., the MG3 paper) show that even at 14 km scale, the inclusion of rimed ice changes the timing and location of precipitation in the western United States due to the different fall speeds and lifetimes of graupel, which is formed when higher vertical velocities result. This effect is expected to be larger at 3 km.

All simulations have been conducted with 58 vertical levels up to 43 km. Set A also includes experiments using 32 vertical levels. We have used the default radiation time step (1 h). The physics and dynamic time steps are set to default at 1800 s for the ∼1^∘ CAM-FV simulation, and this is the default for CAM6 physics for the nominally 1^∘. For 120 km the MPAS dynamic time step is 900 s and the physics time step is 1800 s. We also use 900 s for the 60 km grid-space experiments, scaling it with reduced mesh spacing. The dynamic time step for the MPAS dycore is 20 s for 60–3 km experiments with the physics time step set to 120 s. Instead of using a 20 s time step for the 60–3 km mesh as scaling would imply, we use a 120 s physics time step for the 60–3 km experiments, in part to reduce computational cost and because other studies have shown acceptable results with this physics time step at comparable mesh spacing (e.g., Zeman et al., 2021). We also recognize that the WUS precipitation as the focus of our study is predominantly orographically forced, whereas the physics-time-step-critical processes are related to unstable deep convection, perhaps lending support for a longer physics time step in this application. We acknowledge the possible sensitivity of our results to the physics time step, and we will be examining this more in future work. The average cost for 60–3 km simulations including writes and restarts is ∼4000 to 6000 core hours for 1 d simulation (i.e., ∼120 000 to 180 000 core hours for obtaining 30 d output) using the Cheyenne supercomputer with the scaling of the high-performance computing to be further improved. We would like to acknowledge that model tuning is not performed. Given the interannual variability in precipitation over the WUS study region, we also acknowledge that it is not our goal to reproduce the recent historical climatology but to evaluate the overall model performance.

Table 1 A list of experiments in this study and the key configuration information.

Figure 1 SIMA-MPAS mesh configuration for the 60–3 km experiments. (a) The global-domain mesh configuration with total grid columns of 835 586; (b) the zoomed-in region (see the red box depicted in panel a) for the mesh structure from 60 to 3 km.

2.2 Observations and observationally based gridded products used to evaluate model performance

In this work, we have employed observations from Clouds and the Earth's Radiant Energy System (CERES) Energy Balance and Filled (EBAF) products (Kato et al., 2018; Loeb et al., 2018) for cloud and radiation fluxes properties. We have used GHCN V2 (Global Historical Climatology Network Version 2) gridded data (Fan and van den Dool, 2008) for the land 2 m air temperature globally, which are provided by the NOAA/OAR/ESRL PSL. We have also used Parameter-elevation Regressions on Independent Slopes Model (PRISM) data for gridded observed precipitation and temperature features (Daly et al., 2017) and gridded 4 km observational data for the snow water equivalent (Zeng et al., 2018). We have also used the recently released Livneh precipitation data (Pierce et al., 2021), comprising another gridded observationally based precipitation dataset to better account for extreme precipitation. Another important dataset used for comparison is the WRF (Weather Research and Forecasting) model 4 km simulation data over CONUS from Rasmussen et al. (2021, https://rda.ucar.edu/datasets/ds612.5, last access: January 2022), which used the mean of the CMIP5 model as the boundary forcing. We extracted the same historical time data as the 60–3 km simulations for direct evaluation (i.e., nonhydrostatic CESM vs. nonhydrostatic WRF as a widely used regional climate model).

Detailed descriptions of the open shared datasets used in this study are given below.

• CERES EBAF data products. We use gridded data from the EBAF product from the NASA CERES project, described by Loeb et al. (2018). CERES provides high-quality top-of-the-atmosphere radiative fluxes and cloud radiative effects, as well as consistent ancillary products for the liquid water path (LWP) and cloud fraction. We start with monthly mean gridded products at 1^∘ and make a 20-year climatology from 2000–2020.

• GHCN_CAMS Gridded 2 m Temperature (Land). Global analysis monthly data from the NOAA PSL come with a resolution at 0.5×0.5^∘. They combine two large networks of station observations including the GHCN V2 and the CAMS (Climate Anomaly Monitoring System), together with some unique interpolation methods (https://psl.noaa.gov, last access: January 2022; Fan and Van, 2008).

• PRISM observed data. The PRISM gridded observed data for daily precipitation and daily 2 m air temperature are used at 4 km grid resolution (Daly et al., 2017; https://prism.oregonstate.edu/, last access: January 2022). Covering the continental USA, PRISM takes the station observations from the Global Historical Climatology Network – Daily (GHCN-D) dataset (Menne et al., 2012) and applies a weighted regression scheme that accounts for multiple factors affecting the local climatology (Daly et al., 2017).

• Livneh gridded observationally based precipitation dataset. In addition to PRISM data, to better account for extreme precipitation, recently released Livneh precipitation data (Pierce et al., 2021; http://cirrus.ucsd.edu/~pierce/nonsplit_precip/, last access: September 2022) are also used for model evaluation. The data (∼6 km grid resolution) are shown to perform significantly better in reproducing extreme precipitation metrics than previous version (Pierce et al., 2021).

• Snow water equivalent (SWE) data over CONUS. This is the observational data product we use for snowpack diagnostics. The data are available from the National Snow and Ice Data Center (NSIDC) (at https://nsidc.org/data/nsidc-0719/versions/1, last access: September 2022). The product provides daily 4 km SWE from 1981 to 2021, developed at the University of Arizona. The data assimilated in situ snow measurements from the Snow Telemetry (SNOTEL) network and the Cooperative Observer Program (COOP) network with modeled, gridded temperature and precipitation data from PRISM (Zeng et al., 2018; Broxton et al., 2019).

• CONUS (Continental U.S.) II High Resolution Present and Future Climate Simulation. The WRF (Weather Research and Forecasting) nonhydrostatic model simulations we used for comparison are from Rasmussen et al. (2021) (accessible at https://rda.ucar.edu/datasets/ds612.5, last access: January 2022). The horizontal grid resolution is 4 km with forcing from the mean of the CMIP5 model for both present (1996–2015) and future (2080–2099) mean climate, with hourly output. For the study region we focus on here (i.e., over the western USA), the simulations provide a more realistic depiction of the mesoscale terrain features, critical to the successful simulation of mountainous precipitation (Rasmussen et al., 2021).

The topography details are shown in Fig. 2 over the western US study region, showing that the complex terrains over coastal and mountainous regions have been well resolved in SIMA-MPAS at 3 km resolution (in contrast to 60 km). This is comparable to the topography details in the WRF mesoscale model at a similar resolution, although we do notice the smoother topography in SIMA-MPAS over the 3 km mesh bounds and transient domains (see Fig. S1). For future regional refined applications, we would suggest having a reasonably larger domain area than the study region at the finest resolution to accommodate the noise and instability from mesh transition. When applied, we regridded the SIMA-MPAS model data to the same grid resolution as the PRISM observation and WRF reference data (i.e., 4 km). For the regridded method and procedure, first CAM-MPAS data are remapped from unstructured grids to regular rectilinear lat–long grids at 0.03^∘, and then the rectilinear data are regridded to the same grid spacings as the PRISM using the bilinear interpolation. The orographic gravity wave drag scheme in SIMA-MPAS (used in CESM2(CAM6)) uses a “sub-grid” orography to force the scheme. Sub-grid orography is calculated for each grid cell from a standard high-resolution (1 km) digital elevation model. Thus, the sub-grid orography forcing is small at 3 km, is larger at 60 km, and varies with grid cell size. So, the overall drag should be somewhat similar to the scale but partitioned differently between resolved and unresolved scales.

Figure 2 Topography over the western US region. (a) SIMA-MPAS at 3km refinement; (b) SIMA-MPAS uniform 60 km grid mesh; (c) WRF simulations at 4 km over CONUS.

3 Results

3.1 Mean climatology diagnostics for CESM with MPAS dycore

As the nonhydrostatic dynamical core is coupled to the CESM framework, we would like to understand the mean climate in SIMA-MPAS and how that compares to a standard hydrostatic core (here, using FV), with the experiments described in Table 1. We evaluate the global context of the new formulation of CESM with a nonhydrostatic dynamical core with both 32 and 58 vertical levels. The 58 layer has higher resolution in the planetary boundary layer (PBL) and in the middle and upper troposphere (about 10 additional levels in the PBL and decreasing vertical grid spacing from 1000 m to ∼500 m near the tropopause). Satellite observations are used for comparison as described in the above Sect. 2.2. Simulation results are averaged over the 5-year output under the present-day climatology (with SST and ice forcings from the mean of the period 1996–2005). That means that simulations are forced with the same climatological monthly mean boundary conditions for sea surface temperature and greenhouse gasses every year to reduce interannual variability.

Figure 3 Zonal mean climatology from 5-year simulations with CESM2 and CAM6 physics using different dynamical cores and vertical levels: (a) liquid water path (LWP), (b) ice water path (IWP), (c) cloud fraction, (d) total precipitation rate, (e) land 2 m air temperature, (f) column drop number, (g) shortwave (SW) cloud radiative effect (CRE), and (h) longwave (LW) CRE. Simulations are the default finite-volume (FV) dynamical core with 32 levels (FV L32: solid blue) and 58 levels (FV L58: dashed blue). Also shown are the MPAS dynamical core with 32 levels (MPAS L32: solid red) and 58 levels (MPAS L58). Observations are shown in green for CERES 20-year climatology (from 2000–2020) for LWP, cloud fraction, SW CRE, and LW CRE and for GHCN_CAMS Gridded 2 m Temperature (Land) from 1990–2010 for (e). Shaded values are 1σ annual standard deviations.

Figure 3 indicates that MPAS simulations have a very similar climate to FV simulations. There are some differences in the tropical ice water path in the Southern Hemisphere tropics and some significant differences in subtropical cloud fraction. The climate differences between 32 and 58 levels are also similar between dynamical cores: decreases in the liquid and ice water path at higher vertical resolution. SIMA-MPAS has slight increases in the cloud fraction and precipitation at higher vertical resolution, while SIMA-FV has little change or slight decreases in the cloud fraction. Land surface temperature is well reproduced when ocean temperatures are fixed with both dynamical cores. The column drop number with CAM-MPAS is lower than CAM-FV but more stable with respect to resolution changes. Subtropical shortwave (SW) cloud radiative effect (CRE) and longwave (LW) CRE have higher magnitudes with CAM-MPAS, consistent with a higher LWP and cloud fraction in these regions, yielding better agreement with the meridional CRE structure. When examining the spatial differences (Figs. S2 and S3), we further found that the differences in the wind over the oceans drive differences in aerosols (sea salt) which alter the aerosol optical depth and droplet concentration. The radiative effects come because of cloud fraction changes: high clouds and specifically the ice water path for the long wave and low clouds and the liquid water path for the short wave. The signal in clouds is stronger at 32 levels (L32; Figs. 3, S2), again probably due to larger differences in the PBL, which is better resolved at 58 levels (L58; Figs. 3 and S3). The microphysics is not as directly related to the cloud fraction, which means interaction with the boundary layer turbulence is important. While these changes are easy to spot, they are not that large and are generally well within some of the tuning which is often done during the model development process.

Analysis of the atmospheric wind and temperature structure (Figs. S4 and S5) indicates that SIMA-MPAS compares as well as or better than SIMA-FV with regards to reanalysis winds and thermal structure in the vertical, though biases are different and of a different sign in many regions of the middle atmosphere. There are differences in low-level wind speed and the subtropical jets between MPAS and FV (Fig. S4), driving differences in temperature between them (Fig. S5), particularly in the stratosphere and near the South Pole. The stratosphere and free-troposphere wind differences are due to slightly different damping and deposition of gravity wave drag forcing. The temperature changes above the surface respond to those wind changes. The near-surface temperature differences (e.g., around Antarctica) also relate to transport of air around topography which is different between MPAS and FV.

Overall, SIMA-MPAS produces a reasonable climate simulation, with biases relative to observations that are of similar magnitude to SIMA-FV simulations, despite limited adjustments being made to momentum forcing. SIMA-MPAS has a realistic zonal wind structure with subtropical tropospheric and polar stratospheric jets. There are differences in magnitude from ERA-Interim (Fig S4), but MPAS (which has not been fully tuned) produces a realistic wind distribution. Further tuning of momentum in the dynamical core and physics could reduce these biases. The key feature of this work is that biases in the Northern Hemisphere mid-latitude tropospheric winds are very small for both FV and MPAS. For the temperature profile, there are patterns of bias between the high and low latitudes indicating different stratospheric circulations between the model and the reanalysis. That could be adjusted with the drag and momentum forcing in the model. Note that no adjustment of the physics has been performed.

3.2 Precipitation distribution and statistics

3.2.1 Mean precipitation features

In the western USA during the wet seasons, most of the precipitation occurs over the mountainous regions, with significant impacts on both water resources and potential flood risk management (Hamlet and Lettenmaier, 2007; Dettinger et al., 2011; Huang et al., 2020a). In Fig. 4, we show the wet-season mean (mid-November to mid-March as investigated here) precipitation features over the targeted region with differences from observations. Although the observational differences between PRISM and Livneh on average are small, it provides a more robust evaluation for both mean and extreme precipitation by having those two observational products. The result demonstrates that SIMA-MPAS can well simulate the precipitation intensity and spatial distributions, as compared to PRISM and Livneh observations. The spatial features at 3 km are well captured with a spatial correlation of about 0.93, with precipitation mainly distributed over the Cascade Range, WUS coastal range, Sierra Nevada, and Rocky Mountains. If looking at the precipitation at the coarser resolution (60 km, Fig. S6a) in SIMA-MPAS, the mean domain average of the precipitation (∼2.43 mm, when averaged over years 2000–2002) is similar to the fine-resolution results (∼2.61 mm) but lacking important regional variability and spatial details.

In terms of biases when compared to PRISM data, SIMA-MPAS 3 km overall underestimates the precipitation by about 0.07 mm (bias averaged over the plotted domain), especially over the windward regions, which could relate to the bias in heavy-precipitation frequency and/or the discrepancies in atmospheric river (AR) landfall locations and magnitude from what was observed over the 5-year (wet-season) simulation statistics. We acknowledge that the interannual variability and the sample size of the ARs could also affect the results of landfall precipitation. WRF, on the other hand, tends to overestimate the precipitation in most regions (for about 0.53 mm, bias averaged over the plotted domain compared to PRISM) except for the northwest coast and some Rocky Mountain regions, which can be seen from the relative-difference plot (Fig. 4c). The relative differences in precipitation are generally large over the drier regions in SIMA-MPAS. Overall, compared to PRISM, the bias is negative (for about −0.81 mm on average) over windward regions but positive over the lee side (for about 0.48 mm on average). We also notice that the spatial details of the precipitation are relatively smoothed over the Rocky Mountains, resulting in a large underestimation bias, which could be partly due to the fact that the boundary for the 3 km mesh grids is nearing those regions (see Figs. 1, 2, and S1).

Figure 4 Mean simulated precipitation and differences from observation: (a) wet-season (mid-November to mid-March) daily precipitation intensity over the western USA (1999–2004); (b) absolute differences from PRISM reference; (c) similar to (b) but for relative differences from PRISM (grid box values of less than 1 mm d⁻¹ have been masked) with the SIMA-MPAS model data regridded to the same resolution as the PRISM grid spacings (i.e., 4 km).

Over the western USA, especially in the coastal states, heavy precipitation can be induced by extreme storm events, mainly in the form of atmospheric rivers (Leung and Qian, 2009; Neiman et al., 2011; Rutz et al., 2014; Ralph et al., 2019; Huang et al., 2020b). The capability to capture and predict such extreme events is a significant part of the application of weather and climate models (Meehl et al., 2000; Sillmann et al., 2017; Bellprat et al., 2019). To figure out the performance of SIMA-MPAS in reproducing the precipitation frequency distribution, we combine all the daily data from all the grid points at each coastal state (California, Oregon, and Washington) to calculate the frequency of daily precipitation by intensity (Fig. 5). SIMA-MPAS captures a reasonable distribution of precipitation intensity with respect to PRISM and Livneh observations, with smaller biases than WRF over California and Oregon regions, particularly at more extreme values (such as when daily intensity exceeds 20 mm d⁻¹). We also notice that over the Washington region, the biases for SIMA-MPAS and WRF are at similar magnitudes compared to the observations, although the two observations also show some uncertainties at the upper tail distributions.

Further, when examining the precipitation days with intensity less than 10 to 15 mm d⁻¹, SIMA-MPAS shows a close match to observations, while WRF tends to slightly underestimate the probability. For more extreme precipitation days, models tend to diverge in terms of the behaviors, with SIMA-MPAS showing some underestimation over California and Washington regions (for an average of ∼14 %, ∼7 %, and ∼18 % bias for days when intensity exceeds 20 mm d⁻¹ and is less than 60 mm d⁻¹ for California, Oregon, and Washington respectively). WRF generally overestimates the heavy-precipitation frequency to a much larger extent (for an average bias of ∼42 %, ∼51 %, and ∼18 % for California, Oregon, and Washington respectively). The sign of the biases is consistent with the previously discussed mean precipitation biases. It is not known to us why the biases in SIMA-MPAS are smaller than in WRF. One hypothesis that would limit precipitation intensity is that SIMA-MPAS has strict conservation limits for energy and mass throughout the model, which are not present in WRF. This is a subject for future work but may also be dependent on the specific WRF physics options used. We acknowledge that the initialization without nudging conditions in SIMA-MPAS simulations does not necessarily reproduce monthly or higher time variability but is able to obtain the seasonal means and distributions. We also acknowledge that the interannual variability and the sample size of the ARs could also affect the results of landfall precipitation. Still, those analyses further testify the capability of using SIMA-MPAS for precipitation studies, giving us good confidence in using SIMA-MPAS for storm event studies.

Figure 5 Probability distribution of daily precipitation intensity. All the daily datasets from the five wet seasons for all grid points in each state are used to construct the distribution statistics. The blue lines refer to WRF reference data; the black lines are for the PRISM observation; the dark golden line refers to the Livneh observation; and the SIMA-MPAS results are in red-colored lines. The SIMA-MPAS model data are regridded to the same resolution as the PRISM grid spacings (i.e., 4 km). The x axis starts from 1 mm d⁻¹, and the y axis is transformed with a logarithmic scaling for better visualization of the upper tail distribution.

3.2.2 MG2 vs. MG3 microphysics for simulated precipitation in SIMA-MPAS

We would like to point out that we have used the default microphysics scheme – MG2 (Gettelman et al., 2015) – when configuring those experiments from CESM2. We acknowledge that MG3 (including rimed ice, graupel in this case) could be a better option with the rimed hydrometeors added (see Gettelman et al., 2019), especially when pushing to mesoscale simulations and for orographic precipitation. In detail, Gettelman et al. (2019) found that the addition of rimed ice improved the simulation of precipitation in CESM at 14 km resolution with wintertime orographic precipitation, due to altering the timing of precipitation by more correctly representing the pathways for precipitation formation with more highly resolved scale vertical velocities. To fulfill this caveat but still make the best use of current simulation data, we have conducted another three experiments using the MG3 microphysics scheme for three wet seasons (1999–2002). Similar diagnostics have been performed to those in the previous part but for the results from these three wet seasons (as shown in Fig. 6).

Overall, the precipitation statistics are well represented in SIMA-MAPS compared to observations with both MG2 and MG3 when evaluating from the same three wet seasons. Although still outperforming WRF output, we do recognize that MG2 tends to underestimate heavy-precipitation frequency in certain regions compared to observations, while MG3 produces more intense precipitation with some overestimations over heavily precipitated regions, mostly over the Cascade Range and WUS coastal range (Fig. 6a). From the frequency distributions (Fig. 6b), MG2 and MG3 microphysics both perform well over the study region. Specifically, MG3 produced stronger precipitation than the MG2 output over the Washington region, showing a closer match to the observations than MG2 results. Due to interannual variability, we still need to investigate more different cases, and it is our next-step plan to further investigate the model performance with more test beds.

Figure 6 MG2 vs. MG3 microphysics used in SIMA-MPAS for the wet-season precipitation over the western USA (1999–2002). (a) Mean precipitation intensity; (b) probability distribution of daily precipitation frequency, like Fig. 5 but for three wet seasons with SIMA-MPAS (MG3) added in dashed red lines; again, the SIMA-MPAS model data are regridded to the same resolution as the PRISM grid spacings (i.e., 4 km).

3.3 Accumulated snowpack features

Snowpack characteristics have remained poorly represented in global climate models, lacking high-resolution terrain realization, fine-scale land–atmosphere coupled processes, and interactions with snow's complicated thermal and hydrological properties (DeWalle and Rango, 2008; Liu et al., 2017; Kapnick et al., 2018). Facing this long-standing issue, we expect that with much improved precipitation features and temperature and substantially better resolved complex terrains, snowpack features can be much better represented in CESM. Here, we have compared the accumulated snow water equivalent (SWE) results, which refer to the total accumulated snow from mid-November to mid-March (based on daily output), and then averaged them over the five seasons (see Fig. 7). By comparing these results with the gridded snow water equivalent observational data, it is shown that SIMA-MPAS (MG2) can produce much improved estimation of the snowpack over the mountainous regions, with less overestimation than WRF simulations at similar resolution. However, the overestimation is notable for both SIMA-MPAS and WRF simulations, highlighting the further need to investigate the land–air interactions in rain and snow processes and partitions from the precipitation contribution. In general, SIMA-MPAS can simulate reasonable spatial details for snowpack distribution over mountainous regions (mainly over the Cascade Range, WUS coastal range, Sierra Nevada, and Rocky Mountains) with positive bias over the northern Cascade Range and certain Sierra Nevada mountainous regions.

Figure 7 Wet-season snow water equivalent (SWE) over the western USA. First row: seasonal mean SWE averaged over (1999–2004) from (a) SIMA-MPAS, (b) gridded observation for SWE as described in the Sect. 2.2, and (c) WRF data; Second row (d–f): absolute differences from observation with all data regridded to 4 km for SIMA-MPAS and WRF averaged over 1999–2004 and SIMA-MPAS (MG3) averaged over 1999–2002.

As the snowfall is dominated by the near-surface temperature and precipitation values, we have examined the 2m temperature (T2) here to see how well temperature is captured in SIMA-MAPS. In Fig. 8, the mean T2 (T2mean) is shown averaged over all simulated wet seasons. In general, near-surface temperature results from SIMA-MPAS are overall matched with observations across varied climate zones including coastal, agriculture, desert, inland, and mountainous. However, we also notice that SIMA-MPAS tends to be warmer over most places (with the averaged bias of about 0.65 ^∘C over the plotted domain), except over very high mountain top ranges with cooler bias. On average, the difference for the regions with warmer biases is about 1.35 ^∘C and the difference for those areas with cooler biases is about −0.99 ^∘C when compared to PRISM data. On the contrary, WRF tends to be cooler in most regions except the southern part of Central Valley and some desert regions in the southwestern USA (the average bias is about −1.84 ^∘C over the plotted domain). We have also investigated the T2 bias in the 120 km simulations to see if this is a consistent model bias. By comparing FV and MPAS together (Fig. S7), it turns out that SIMA-MPAS tends to be warmer with higher net surface shortwave and longwave fluxes over the wet-season period discussed here (Fig. S8). Still, overall, the land model coupled with the atmosphere also does a good job here under a realistic topography.

Figure 8 Daily mean 2 m air temperature (T2mean) averaged over 1999–2004, November–March. (a) PRISM observation dataset; (b, c) the differences between SIMA-MPAS and PRISM and between WRF and PRISM respectively (note for the difference plots, all data are regridded to the same resolution as PRISM).

3.4 Large-scale moisture flux and dynamics

Further, we have investigated the wind profile that directly connects to the subtropical to middle-latitude moisture fluxes over the northeast Pacific and hits western US regions. First, we have examined the cross-sections of zonal and meridional wind patterns (at 130^∘ W, near the western US coast) at both 60–3 and 60 km to determine the dynamic changes with the refinement mesh (Fig. 9). As we can see, the mean westerly zonal winds are about 10 % stronger at the jet stream level near 200–250 hPa in 60–3 km simulations compared to the 60 km results. The mean meridional wind (dominantly southward) however is weaker in 60–3 km simulations than in the 60 km ones. The precipitation over the western US coast is largely associated with the concentrated water vapor transport over the North Pacific, known mainly in the form of atmospheric rivers (Rutz et al., 2014). It is our intention to investigate the wind dynamics transitioning from the coarse scale to mesoscale in future work. As another source of the precipitation uncertainty, we would like to acknowledge the sensitivity from the physics time step (see Fig. S9) when comparing the precipitation in 60–3 km simulations (a shorter physics time step) to the 60 km results at the regions with the same grid resolutions.

Figure 9 Composite wind profile along the western US coast (cross-section at 130^∘ W, near the western US coast) (averaged over 2000–2002, November–March). (a, c) Mean latitude-height cross-section of zonal winds (m s⁻¹) for SIMA-MPAS 60–3 km (a) and 60 km (c); (b, d) similar to (a) and (c), except for meridional winds.

In Fig. 10, we further examine the large-scale moisture flux pattern from the integrated water vapor transport (IVT) in the set of simulations with and without regional refinement. The spatial pattern of the moisture flux is generally similar between those two sets of experiments, dominated by the zonal winds (see Fig. 9). If checking the IVT values along the longitude of 130^∘ W, the differences (about 3 % on average) are quite small along the WUS extent. With the large-scale dynamics and local fine-scale processes well integrated into this nonhydrostatic global climate model, it gives confidence in reproducing and predicting precipitation across the weather and climate scales.

Figure 10 Mean instantaneous vertically integrated water vapor flux transport over the western USA (2000–2002, November–March): (a) SIMA-MPAS 60–3 km and (b) SIMA-MPAS 60 km. Wind is overlaid for the averaged lower levels (height from ∼500 to ∼2000 m).

4 Summary and discussion

In this study, we describe SIMA-MPAS, which is built upon the open-source Community Earth System Model (CESM) with a nonhydrostatic dynamical core, the Model for Prediction Across Scales (MPAS). We would like to try to answer several questions about the performance of this new-generation model when applying at convection-permitting resolutions and when bridging both weather- and climate-scale simulations in a single global model. We have chosen the western USA as our study region to examine the precipitation features in SIMA-MPAS at fine scales and how the model performs when compared to both observations and a regional climate model.

To answer those questions, we have designed and conducted a set of experiments. First, we have tested CESM at the same coarse resolution using both MPAS as the nonhydrostatic core and finite volume as the hydrostatic core for multiple years of climatology. Secondly and as the focus of this work, a variable-resolution mesh is configured with 3 km refinement centered over the western USA. We have performed five separate wet-season simulations to obtain the precipitation statistics. In addition, we have also included uniform 60 km simulations from the model for two seasons.

We first evaluated the mean climate in SIMA-MPAS to see how that compares to the hydrostatic model counterpart (here, SIMA-FV). The diagnostics show that MPAS simulations have a very similar climate to FV simulations. SIMA-MPAS has slight increases in cloud fraction and precipitation at the higher vertical resolution, while SIMA-FV has little change or slight decreases in cloud fraction. Overall, SIMA-MPAS produces a reasonable climate simulation, with biases relative to observations that are not that different from SIMA-FV simulations, despite limited adjustments being made to momentum forcing, and no adjustment of the physics has been performed.

When compared to both observations and a traditional regional climate model at similar fine resolutions for mean precipitation and heavy-precipitation behaviors, SIMA-MPAS can capture the spatial pattern and mean intensity (with the spatial correlation of about 0.93 relative to PRISM), which is also comparable to WRF results. We do notice there are some underestimations, mostly in SIMA-MPAS, and overestimations, mostly in WRF. Further, SIMA-MPAS captures the distribution of precipitation intensity with respect to observations with smaller biases than WRF over California and Oregon regions, particularly at more extreme values. With additional experiments, SIMA-MPAS with MG3 microphysics (graupel) produces stronger precipitation than the MG2 version (as used in other experiments in this study as the default microphysics scheme), and the MG3 results also well represented the precipitation statistics for both spatial mean and frequency distribution. The difference between MG3 and MG2 is the rimed hydrometeors added to MG3 (see Gettelman et al., 2019, for detailed descriptions), which could matter more when pushing to mesoscale simulations and for orographic precipitation. We also acknowledge the interannual variability, and it is our next-step plan to further investigate the model performance with more test beds.

We further show that SIMA-MPAS can produce much improved estimation of the snowpack over the mountainous regions compared to coarse resolutions, with less overestimation than WRF simulations at similar resolution. In general, SIMA-MPAS can simulate some reasonable spatial details for snowpack distribution over mountainous regions (mainly over the Cascade Range, WUS coastal range, Sierra Nevada, and Rocky Mountains) with positive bias over the northern Cascade Range and certain Sierra Nevada mountainous regions. The overestimation is notable for both SIMA-MPAS and WRF simulations, needing further investigations. We also notice that SIMA-MPAS tends to be warmer over most places, except over very high mountain top ranges with cooler bias.

The results further testify to the capability of using SIMA-MPAS for precipitation studies, giving us good confidence in using SIMA-MPAS for storm event studies. We focus on multiple-season statistics for model performance. Given the large-scale dynamics and local fine-scale processes well integrated into this nonhydrostatic global climate model, it shows promise in reproducing and predicting precipitation across the weather and climate scales. It is our further intention to investigate the wind dynamics transitioning from the coarse scale to mesoscale in future work and to further investigate the model performance with more test beds for convection-permitting weather and climate systems across scales.