Evaluation of a unique approach to high-resolution climate modeling using the Model for Prediction Across Scales – Atmosphere (MPAS-A) version 5.1

We present multi-seasonal simulations representative of present-day and future environments using the global Model for Prediction Across Scales – Atmosphere (MPAS-A) version 5.1 with high resolution (15 km) throughout the Northern Hemisphere. We select 10 simulation years with varying phases of El Niño–Southern Oscillation (ENSO) and integrate each for 14.5 months. We use analyzed sea surface temperature (SST) patterns for present-day simulations. For the future climate simulations, we alter present-day SSTs by applying monthly-averaged temperature changes derived from a 20-member ensemble of Coupled Model Intercomparison Project phase 5 (CMIP5) general circulation models (GCMs) following the Representative Concentration Pathway (RCP) 8.5 emissions scenario. Daily sea ice fields, obtained from the monthly-averaged CMIP5 ensemble mean sea ice, are used for present-day and future simulations. The present-day simulations provide a reasonable reproduction of large-scale atmospheric features in the Northern Hemisphere such as the wintertime midlatitude storm tracks, upper-tropospheric jets, and maritime sea-level pressure features as well as annual precipitation patterns across the tropics. The simulations also adequately represent tropical cyclone (TC) characteristics such as strength, spatial distribution, and seasonal cycles for most Northern Hemisphere basins. These results demonstrate the applicability of these model simulations for future studies examining climate change effects on various Northern Hemisphere phenomena, and, more generally, the utility of MPAS-A for studying climate change at spatial scales generally unachievable in GCMs.

1 Introduction

We present a novel approach to high-resolution climate modeling with the intent of examining the effects of climate change on high-impact Northern Hemisphere weather phenomena. It is nearly certain that rising global greenhouse gas concentrations over the next century will result in significant changes to the Earth's climate system (IPCC, 2014). Further understanding of how climate change will affect global and regional weather is essential to informing the scientific community, stakeholders, and policymakers on what actions should be taken to prepare for the future. Here, we present a unique set of high-resolution, multi-seasonal, global atmosphere-only simulations conducted with the Model for Prediction Across Scales – Atmosphere (MPAS-A; Skamarock et al., 2012) in present and future environments for the purpose of studying climate change effects on Northern Hemisphere weather phenomena, including extreme events. Present-day conditions are simulated using the current atmospheric composition and observed lower boundary conditions; climate change is represented by modifying the atmospheric composition to values appropriate for the late 21st century, and applying consistent changes in lower boundary conditions derived from CMIP5 coupled global climate models. Through its variable-resolution grids, MPAS-A offers the possibility of investigating local weather phenomena at high resolution in the context of a global model, while avoiding the prohibitive demands on computational resources entailed by running a model globally at high resolution. To our knowledge, however, climate change experiments at long integration times are a novel application of MPAS-A. Therefore, in order to demonstrate their utility for addressing climate change effects on high-impact weather events, it is necessary to evaluate how large-scale circulations and responses to warming are represented in such simulations, thus defining the objective of this paper.

With simulations spanning several centuries and multiple ensemble members, and the inclusion of atmosphere–ocean coupling, the latest generation of general circulation models (GCMs) from the Coupled Model Intercomparison Project phase 5 (CMIP5) are common tools for determining the effects of climate change. Due to current computational limitations, however, the grid spacing of these simulations is largely restricted to ∼1^∘ (∼100 km) or greater. While this coarse resolution is suitable for representing large-scale atmospheric features such as the polar amplification of global warming and teleconnections, it is insufficient for resolving weather extremes, especially those associated with smaller-scale systems such as tropical cyclones, mesoscale features within extratropical cyclones, and convective storms (e.g., Mizielinski et al., 2014 and references therein; Small et al., 2014; Prein et al., 2015; Haarsma et al., 2016; Roberts et al., 2018). These smaller-scale systems often result in significant socioeconomic impacts; therefore, in order to fully ascertain the societal impacts of climate change, it is essential to complement existing GCM simulations with simulations at resolutions sufficiently fine to capture these high-impact phenomena. The ongoing High Resolution Model Intercomparison Project (HighResMIP; Haarsma et al., 2016) associated with the CMIP6 will also be highly beneficial in understanding how weather extremes respond to climate change.

To date, several model downscaling techniques have been developed for this purpose. For example, regional downscaling (e.g., Wang et al., 2004; Giorgi et al., 2009) computationally allows for finer grid spacings by employing a smaller domain, thus circumventing the resolution deficiency of traditional GCMs. Using a regional domain, however, presents the issue of how to specify lateral boundary conditions, and two-way interactions with larger scales cannot be fully incorporated (Small et al., 2014). Global models eliminate the constraints of lateral boundaries but are expensive to run for long periods at high resolutions. Incorporating nests within a global domain, or using mesh refinement grids, however, can help alleviate this expense.

Another useful method for assessing climate change effects is the “pseudo-global warming” (PGW) method, initially called “surrogate global warming” (e.g., Schär et al., 1996; Frei et al., 1998; Kimura and Kitoh, 2007; Hara et al., 2008; Rasmussen et al., 2011; Mallard et al., 2013; Lackmann, 2013, 2015; Trapp and Hoogewind, 2016). In PGW experiments, high-resolution control simulations are conducted, typically replicating an observed weather event. The high-resolution initial and boundary conditions are then modified with “delta” fields derived from GCMs, and the event is re-simulated, allowing assessment of changes in the characteristics of the event as a function of larger-scale environmental change. An important advantage of the PGW method is that realistic, high-resolution synoptic-scale and mesoscale settings are guaranteed. This method is consistent with the “storyline” approach described by Shepherd (2016), Hazeleger et al. (2015), and Trenberth et al. (2015). A limitation of PGW case studies is the inability to study the frequency of occurrence of such events. To alleviate this limitation, some investigators have conducted long-duration regional PGW simulations (e.g., Ban et al., 2014; Willison et al., 2015; Liu et al., 2017), which allow for analysis of statistical changes extending beyond the case study of a single event. All regional PGW experiments, however, are limited by the need to impose lateral boundary conditions, which reduces the dynamical freedom of the simulations.

Given recent advances in computational power and data storage, several modeling groups have performed long-term high-resolution global simulations, both with atmosphere-only and coupled atmosphere–ocean configurations (e.g., Small et al., 2014; Kodama et al., 2015; Murakami et al., 2015; Roberts et al., 2015 and references therein; Haarsma et al., 2016; Yamada et al., 2017). Models that include coupling between the atmosphere and ocean have the advantage of two-way communication, allowing the possibility of realistic atmosphere–ocean interactions. At long integration times, however, climatologies of coupled models have been known to suffer from biases due to the drift in sea surface temperatures (SSTs), which can negatively affect regional climate projections (e.g., He and Soden, 2016). Previous studies have determined that resolution is an important factor for a more accurate representation of synoptic and mesoscale phenomena in the atmosphere (Willison et al., 2013; Small et al., 2014; Prein et al., 2015) and thus should be maximized whenever possible.

Figure 1 Variable-resolution mesh for MPAS-A simulations and geographical regions of the tropical cyclone basins defined in this study.

The present global simulations use a 15 km grid in the Northern Hemisphere, relaxing to a 60 km grid in the Southern Hemisphere to reduce computational expense (Fig. 1). These simulations are conducted with MPAS-A as individual time-slice runs, selected to span a range of El Niño–Southern Oscillation (ENSO) states. We use high-resolution SST analyses that capture oceanic eddies and fronts, which have been shown to exert an important influence on atmospheric variability (e.g., Kirtman et al., 2012; Siqueira and Kirtman, 2016; Ma et al., 2017; Parfitt et al., 2017). Our novel modeling approach aims to eliminate several limitations from the previously discussed methods. For one, using a global model circumvents issues related to lateral boundary conditions, thus improving upon limited-area simulations; here, higher resolution in the Northern Hemisphere was obtained at the expense of reduced resolution in the Southern Hemisphere. Our future simulations are similar to PGW in the treatment of SST in that we apply a GCM-based delta to analyzed SST fields; this incorporation of high-resolution SSTs precludes issues of the type noted by He and Soden (2016), and represents a potential improvement to coupled atmosphere–ocean model configurations. The atmospheric resolution of 15 km is sufficiently high to represent strong tropical and extratropical cyclones as well as flooding rainfall, which is a considerable improvement compared to GCM simulations. Furthermore, the sample size is sufficient to allow statistical comparisons of features such as Northern Hemisphere storm tracks and tropical cyclone activity, thus improving upon the traditional PGW case study approach.

We present our simulations with the intention of providing an additional realization of a complex system in order to improve our understanding of potential climate change effects on Northern Hemisphere high-impact weather phenomena. For such a modeling system to be useful for this purpose, it is necessary that

• it reproduces the present-day climate and global circulation of the atmosphere;

• it demonstrates the benefits of enhanced resolution in simulating high-impact weather phenomena; and

• it provides simulations of a future climate consistent with expectations derived from GCMs.

Thus, in the present paper, following a discussion of the model and how our simulations were conducted (Sect. 2), we offer analyses of its representation of the present-day climate (Sect. 3), and its simulation of Northern Hemisphere tropical cyclones and their climatology within the present-day climate (Sect. 4). We focus on tropical cyclones as an exemplar of high-impact weather phenomena that are challenging to represent accurately in models, and for which successful simulation demands high resolution (e.g., Davis, 2018). In Sect. 5, we examine the model representation of climate change in response to global warming boundary conditions. Last, Sect. 6 presents a summary of our findings and discusses future applications of our simulations for investigating how high-impact weather may change in a warmer climate. These applications include ongoing research efforts investigating climate change effects on the extratropical transition of tropical cyclones (TCs), TC seasonality, and persistent anomalies and blocking, all of which will be subjects of future publications.

2 Models, experiments, and performance

2.1 Model configuration

We conduct our simulations using the atmospheric component of MPAS-A, version 5.1 (Skamarock et al., 2012). MPAS-A is a non-hydrostatic global atmosphere-only model that uses unstructured Voronoi meshes (Du et al., 1999) to create variable-resolution grids. This grid structure permits localized areas of high resolution to transition gradually to lower resolutions, thus alleviating the boundary issues associated with sharp transitions between domains in traditional nesting approaches (Park et al., 2014). The focus of the simulations presented here is on Northern Hemisphere phenomena; therefore, we use a variable-resolution mesh with 15 km grid spacing over the Northern Hemisphere, expanding out to 60 km in the Southern Hemisphere (Fig. 1).

The MPAS-A atmospheric physics suite includes a subset of schemes adapted from versions of the Weather Research and Forecasting (WRF) model (Skamarock et al., 2008). Our simulations employ the following physics parameterizations: WRF Single-Moment 6-Class Microphysics (WSM6; as in WRF 3.8.1), Yonsei University (YSU; as in WRF 3.8.1) representation of the planetary boundary layer, Tiedtke (as in WRF 3.3.1) subgrid-scale convective parameterization, Community Atmosphere Model (CAM; as in WRF 3.3.1) shortwave and longwave radiation, and the Noah land surface model (as in WRF 3.3.1) for surface processes. We selected the Tiedtke convective parameterization scheme because it includes convective momentum transport (CMT), which has been shown to be important for reducing model biases in surface winds and TC intensity (Zhang and McFarlane, 1995; Han and Pan, 2006; Hogan and Pauley, 2007; Richter and Rasch, 2008). CMT also improves the representation of features such as the Intertropical Convergence Zone (Zhang and Wang, 2006; Kim et al., 2008). We completed a series of preliminary tests using a quasi-uniform 60 km mesh to further refine our physics choices (not shown).

2.2 Present-day and future climate simulations

We selected 10 simulation years with varying phases of ENSO based on the Multivariate ENSO Index (MEI) and the Oceanic Niño Index (ONI) over the TC season (Table 1). These years were also chosen to sample a range of TC activity in the North Atlantic, eastern North Pacific, and western North Pacific basins. We chose to sample phases of ENSO rather than other modes of climate variability due to its strong connection to global TC activity (e.g., Gray, 1984; Chan, 1985; Lander, 1994; Chu and Wang, 1997; Kossin et al., 2010). Each simulation is integrated for 14.5 months, from 1 March of the first year through 14 May of the following year, with the first month of each simulation discarded as spin-up; output is recorded every 6 h.

Table 1 Average Multivariate ENSO Index (MEI), Oceanic Niño Index (ONI), and corresponding ENSO phase during the TC season (June–November) for the chosen simulation years.

Table 2 List of 20 CMIP5 GCMs used to compute ensemble mean temperature “deltas” and sea ice fields.

We used the ECMWF ERA-Interim Reanalysis (ERA-I; European Centre for Medium-Range Weather Forecasts, 2009; Dee et al., 2011) with a spectral T255 resolution (∼0.7^∘ horizontal grid spacing) for present-day initial conditions. SST and sea ice fields are updated daily throughout the simulations. The configuration of these surface fields is discussed further in Sect. 2.3. For the future climate simulations, we modify the ERA-I initial and lower boundary conditions by adding monthly-averaged temperature changes derived from a 20-member ensemble of CMIP5 GCMs (Table 2). A similar change was applied to the deep-soil temperature. These temperature changes are calculated by subtracting the 1980–1999 average temperature from the 2080–2099 average temperature following the Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report (AR5) Representative Concentration Pathway (RCP) 8.5 emissions scenario, interpolated to the ERA-I grid and added to the existing temperature data at all atmospheric pressure and soil levels. Geopotential height and specific humidity are adjusted by the model based on the imposed temperature changes; relative humidity is held constant at the initial time in the spin-up run (see the following paragraph). We set carbon dioxide (CO₂) concentrations in the future climate simulations to 936 ppm, the level projected by the RCP8.5 emissions scenario for 2100 (Meinshausen et al., 2011). Present-day CO₂ concentrations are based on analyzed values set according to the year.

Figure 2 Flowchart depicting the “daisy-chain” simulation method.

Rather than running the simulations in chronological order, the simulation years are sorted from the strongest La Niña year (i.e., the year with the smallest MEI and most negative ONI; Table 1) to the strongest El Niño year (i.e., the year with the largest MEI and most positive ONI; Table 1). This design aims to minimize model spin-up in response to changes in SST. With the present-day and future initial conditions set, we conduct full simulations for a neutral ENSO year (e.g., 2013) for both present-day and future environments. These single-year simulations are used as spin-up and are therefore excluded from our analysis. While we took this precaution to allow the model atmosphere to come into equilibrium with the imposed warming and adjusted CO₂ for the future climate experiment, we repeated this process for the present-day simulation to maintain consistency. We then used the output from 1 March, towards the end of the initial spin-up simulation, to initialize the first simulation year. This method continues for both the present-day and future experiments by using the output from the latter part of one simulation (e.g., 1 March) to initialize the next (Fig. 2). Applying this unique “daisy-chain” technique avoids the need for excessive spin-up times for each year; instead, we discard only the output from the first month, which allows any discontinuities arising from the change in SST to equilibrate.

2.3 Lower boundary conditions

As with the atmospheric initial conditions, the SST fields used in the simulations are taken primarily from ERA-I. The SSTs in ERA-I have, however, been derived from several different datasets over the years (Dee et al., 2011). For reanalysis times after February 2009, ERA-I surface fields originate from the Operational Sea Surface Temperature and Sea-Ice Analysis (OSTIA; Donlon et al., 2012). To maintain consistency between all simulations, the OSTIA SST, interpolated from its native 0.05^∘ horizontal grid spacing to the ERA-I grid, is used for simulation years prior to 2009. Therefore, we effectively use OSTIA SST for all simulations. For present-day soil temperature and moisture, we use the ERA-I fields.

Figure 3 Example SST (K) on 1 March 2013 for (a) present-day and (b) future MPAS-A simulations and example sea ice fraction on 1 September for the (c) present-day and (d) future MPAS-A simulations. Contours are shaded every 1 K in panels (a) and (b) and every 0.05 units in panels (c) and (d).

The SSTs for the future climate simulations are altered in the same manner as the initial condition atmospheric and soil temperatures (e.g., Fig. 3a–b). The same technique of adding a GCM delta field onto existing data cannot, however, be used for sea ice. Instead, similar to Mizielinski et al. (2014), monthly-averaged CMIP5 ensemble mean sea ice fields are temporally interpolated to create daily sea ice fields for both present-day (1980–1999) and future (2080–2099 under the RCP8.5 emissions scenario) time periods. An example of these sea ice fields is shown in Fig. 3c–d. We then replaced the analyzed sea ice in the ERA-I with these climatological fields for use in all model simulations. While the climatological present-day sea ice does not entirely match the analyzed field in the ERA-I (e.g., the sea ice edge is much more diffuse), handling the sea ice in this manner ensures that it is plausibly represented in the future climate simulations. The presence of an overly diffuse ice edge could result in unrealistically weak lower tropospheric baroclinicity during warm seasons in these locations.

Our technique for simulating a future climate is similar to the PGW approach in the sense that (1) the analyzed initial and lower boundary conditions are altered by adding projected temperature changes from GCMs to represent future conditions, and (2) analyzed high-resolution SST fields are used to preserve realistic representation of ocean eddies and SST gradients. High-resolution SST is of demonstrated importance for midlatitude cyclone development and other regional climate changes (e.g., Brayshaw et al., 2011; Booth et al., 2012; Kirtman et al., 2012; He and Soden, 2016; Siqueira and Kirtman, 2016). By using a global model, however, one of the main limitations of PGW, the constraint of the lateral boundary conditions, is alleviated. Therefore, our future climate simulations are best described as time-slice experiments with prescribed high-resolution SSTs. The UPSCALE (UK on PRACE: weather-resolving Simulations of Climate for globAL Environmental risk) experiments described by Mizielinski et al. (2014) use a similar time-slice technique for simulating a future climate. By simulating a small ensemble of 26 years, UPSCALE samples a broad range of interannual variability and ENSO states; however, 25 km grid spacing is insufficient for resolving full-strength tropical cyclones (Davis, 2018). Therefore, our simulations complement UPSCALE by offering sufficiently high resolution to better capture the atmospheric mesoscale, specifically tropical cyclones.

2.4 Computational performance

We conducted the MPAS-A simulations on the US National Center for Atmospheric Research (NCAR) supercomputer, Cheyenne (Computational and Information Systems Laboratory, 2017). Cheyenne is a 5.34 petaflop SGI ICE XA cluster with 145 152 Intel Xeon processor cores and 313 terabytes (TB) of memory. Each 14.5-month simulation was run on 1152 cores and consumed roughly 92 000 CPU hours, including resources needed for post-processing, leading to a total of ∼1.9 million core hours used for these experiments.

We post-process model output to vertically interpolate fields to selected isobaric levels and horizontally interpolate from the native unstructured mesh to a $\mathrm{0.15}{}^{\circ }\times \mathrm{0.15}{}^{\circ }$ latitude–longitude grid. Due to storage constraints, we saved a limited number of variables for the Northern Hemisphere only; however, monthly restart files are archived, enabling replication of a particular period of time or event as needed. The post-processed output occupies approximately 50 TB of storage space for the output for all 20 simulations and is currently stored on Cheyenne's High Performance Storage Space (HPSS) and at North Carolina State University.

2.5 Assumptions and limitations

For our future climate simulations, we computed temperature delta values using the mean IPCC AR5 RCP8.5 emissions pathway. While other plausible scenarios exist, we selected a high emission pathway to maximize the signal of climate change in our simulations. Using the GCM ensemble mean temperature changes to alter our initial and lower boundary conditions diminishes the considerable amount of variability in the temperature changes projected by individual GCMs. Computing an ensemble mean from a set of simulations using temperature change fields from each GCM is, however, unlikely to produce significantly different results (Hill, 2010; Lackmann, 2015; Marciano, 2014). An alternate strategy to using the GCM ensemble mean SST change would be to apply changes computed as a function of GCM ENSO phase. However, given uncertainties in the GCM SST change fields, we felt that the use of the ensemble mean was best for our initial set of experiments. Likewise, our analysis does not address possible changes in the distribution of ENSO phases.

The adjustment of geopotential height based on the imposed temperature changes for the future climate simulation introduces some degree of imbalance between the model mass and wind fields. In previous studies, we utilized the digital filter initialization (DFI) capability of the WRF model to reduce these imbalances. Since this feature is not available in MPAS-A, we conducted a full 14.5-month spin-up simulation to allow time for the dynamics of the model atmosphere to restore balance. We also maintain constant relative humidity between the present-day and future simulations in the initial conditions. While this assumption may be appropriate over ocean basins, it does not necessarily hold true over land areas (e.g., Sherwood and Fu, 2014). With no constraints on the lateral boundaries, however, this constraint is only applied once, at the beginning of our spin-up simulation, and therefore is not enforced during the subsequent simulations. The relative humidity within our model domain evolves freely through the duration of the simulations; by the end the initial year-long spin-up period, we expect the distribution of water vapor to be fully equilibrated with the simulated future climate.

While our treatment of sea ice in the model allows for a plausible representation of future conditions, we use identical sea ice fields in each member of our present-day simulation set and similarly for the future set. We therefore exclude the effects of interannual variability in sea ice. Several studies have highlighted the connection between sea ice variability and atmospheric circulations in the Northern Hemisphere (e.g., Deser et al., 2000; Overland and Wang, 2010); our intention here, however, is to minimize this influence and instead focus on changes due to altered temperatures and atmospheric composition. Another limitation inherent in our methods is the assumption that future patterns of SST variability will remain similar to what they are today. Nevertheless, we believe the benefits of using high-resolution SST analyses to preserve realistic SST gradients and alleviate regional biases associated with atmosphere–ocean coupling (e.g., He and Soden, 2016) outweigh this limitation.

Many previous studies have shown that neglect of SST cooling due to cyclone passage results in TCs that are too strong (e.g., Schade and Emanuel, 2009). Use of analyzed SST fields in our simulations does not allow for TC-generated cold wakes, which could contribute towards a positive bias in TC intensities and could lead to unrealistic temporal clustering of TCs. The use of convective parameterization, however, particularly the Tiedtke scheme which adjusts momentum, tends to weaken TCs through momentum adjustment in a warm-core cyclonic structure, an effect opposite to that resulting from the neglect of SST cooling. Additionally, the presence of pre-existing cold wakes in the OSTIA SST field could erroneously weaken TCs that occur in their path; however, these pre-existing cold wakes in the OSTIA do not appear to be particularly strong in magnitude (not shown) and therefore are unlikely have a substantial impact on simulated TC strength. Ideally, a grid length of 4 km or less would be used to fully capture TC structure and intensity (e.g., Gentry and Lackmann, 2010), but computational expense does not allow this for the Northern Hemisphere region of interest for the simulation durations necessary to obtain statistically meaningful results regarding the impacts of climate change. A benefit of our configuration is that the resolution is sufficiently high to capture nearly the full range of TC intensity; preliminary testing highlighted the capability of our 15 km grid to replicate realistic TC structures, including spiral rain bands and a defined eye (not shown). We acknowledge that the neglect of sea-surface cooling, the use of parameterized convection, and potential effects of pre-existing cold wakes in the SST analysis data are limitations to our approach. These limitations are, however, consistent between present-day and future simulations, allowing any differences found in TC intensity to remain meaningful (Patricola and Wehner, 2018).

We recognize that the methods employed in this study account only for projected changes due to increased anthropogenic greenhouse gases and therefore do not represent other external climate forcings. Changes in other aspects of the climate system, such as changes in aerosols, deep soil moisture, and vegetation, are not represented. Despite the limitations discussed, our method alleviates limitations associated with regional PGW and coarse GCMs, and is much more computationally efficient than running high-resolution global models for long integration periods (e.g., centuries); the result is a set of controlled simulations suitable for examining the effects of climate change on high-impact weather events.

3 Model climate: precipitation and midlatitude features

3.1 Extratropical storm tracks

There are two primary midlatitude storm track regions in the Northern Hemisphere (the North Pacific and the North Atlantic) where baroclinic waves form over regions of enhanced temperature contrast linked to warm western boundary currents off the east coasts of Asia and North America and propagate eastward through downstream development (Chang et al., 2002). The extratropical cyclones in these regions play an essential role in the Earth's climate system and contribute to everyday weather, including high-impact events. Therefore, it is important that they are well represented in model simulations.

Figure 4 Average wintertime (DJF) sea-level pressure (SLP) variance (hPa²) over the 10 simulation years for (a) present-day MPAS-A simulations, (b) ERA-I 10-year climatology, and (c) the model bias (MPAS-A minus ERA-I). Contours are shaded every 10 hPa² in panels (a) and (b) and every 5 hPa² in panel (c). MPAS-A output was linearly interpolated to the ERA-I grid for point-to-point comparison. The pattern correlation coefficient is reported in the top right of panel (b). Stippling in panel (c) indicates locations where the MPAS-A 10-year mean exceeds the range computed from 100 random samples of 10-year means from ERA-I by more than 5 % (∼7 hPa²).

As suggested by Chang and Fu (2003), variance in daily-mean fields can be used as proxies for storm track activity. Here, we use the 24 h variance of daily-mean sea-level pressure (SLP), calculated using Eq. (2) from Chang et al. (2013):

$$\begin{array}{}\text{(1)}& \mathrm{SLP}\mathrm{variance}=\stackrel{\mathrm{‾}}{{\left[\mathrm{SLP}\left(t+\mathrm{24}\mathrm{h}\right)-\mathrm{SLP}\left(t\right)\right]}^{\mathrm{2}}},\end{array}$$

where the overbar indicates the quantity is averaged over time, in this case over the winter season (December–February; DJF) when storm activity in the Northern Hemisphere is maximized (Chang et al., 2002; Brayshaw et al., 2009). Figure 4 shows the wintertime SLP variance for the MPAS-A simulations compared to the ERA-I; the ERA-I climatology in Fig. 4b is computed using only the 10 years corresponding to our simulations. The North Pacific and North Atlantic storm track regions are clearly evident in the model simulations; the overall spatial correlation coefficient is greater than 0.98, indicating that general patterns of SLP variance are well reproduced in the MPAS-A simulations. As evident from the positive biases at ∼35^∘ N, 165^∘ W and ∼40^∘ N, 30^∘ W (Fig. 4c), both storm track regions are shifted equatorward, and the North Pacific storm track is more zonally oriented in the MPAS-A simulations. Comparison with 100 random samples of 10-year means from the ERA-I record indicates that these biases, primarily the shift in the North Pacific, likely represent true differences between the MPAS-A simulations and the real atmosphere. Negative biases in simulated storm activity occur east of Greenland, over Scandinavia, and throughout central North America; these differences, however, fall within the range of observed variability (Fig. 4c).

3.2 Northern Hemisphere jet and sea-level pressure features

Corresponding to the Northern Hemisphere extratropical storm tracks are the midlatitude jet features, represented by the wintertime average zonal wind speed in Fig. 5. As in Fig. 4, the ERA-I climatology in Fig. 5b includes only the 10 simulation years. As indicated by a pattern correlation coefficient of ∼0.99, the orientation and spatial extent of the North Pacific and North Atlantic jets at the 250 hPa level are well replicated by MPAS-A (Fig. 5a–b). While the North Atlantic jet maximum is slightly stronger in the MPAS-A simulations, the general strength of both features compares well between the simulations and reanalysis (Fig. 5a–b). Furthermore, examination of a cross-section of zonally averaged zonal wind shows the jet maximizes at roughly the same altitude and latitude (∼200 hPa and ∼30^∘ N) in both the MPAS-A simulations and the ERA-I climatology, albeit the MPAS-A maximum is moderately weaker (Fig. 5c–d). Additionally, semi-permanent maritime SLP features, such as the Aleutian Low over the Bering Sea, the North Pacific subtropical high, and the Icelandic Low, are well captured in the MPAS-A simulations (Fig. 5a–b). The Bermuda High in the North Atlantic, however, while evident in the MPAS-A simulations, is comparatively weaker than analyzed.

Figure 5 (a, b) Average wintertime (DJF) 250 hPa zonal wind speed (m s⁻¹; shaded every 2.5 m s⁻¹) and SLP (hPa; contoured every 4 hPa) over the 10 simulation years for (a) present-day MPAS-A simulations and (b) ERA-I 10-year climatology. (c, d) Average wintertime (DJF) cross-section of zonally averaged zonal wind speed (m s⁻¹; shaded every 2.5 m s⁻¹) for (c) present-day MPAS-A simulations and (d) ERA-I 10-year climatology. The pattern correlation coefficient for 250 hPa zonal wind speed is reported in the top right of panel (b). SLP contours in panels (a) and (b) are masked over land due to noise in areas of complex terrain.

3.3 Large-scale precipitation

Average precipitation over the 10 simulation years compared to the 19-year (1998–2016) climatology from the Tropical Rainfall Measurement Mission (TRMM; Huffman et al., 2007; Tropical Rainfall Measuring Mission, 2011) 3B42 product is shown in Fig. 6. Because 4 of our 10 simulation years occur before the TRMM record began, we opted to use the full TRMM climatology for comparison. A pattern correlation coefficient of ∼0.95 indicates that MPAS-A simulates the general spatial pattern of tropical precipitation well; the Intertropical Convergence Zone (ITCZ) in the equatorial Pacific and maxima along the west coast of India, over the Himalayas, and throughout northern South America are all well-represented by the model. The primary difference between the two precipitation fields is the overproduction of precipitation by MPAS-A in many areas (Fig. 6c), an issue common among other high-resolution modeling studies (e.g., Bacmeister et al., 2014; Small et al., 2014). The overestimation in the subtropical Pacific basin (Fig. 6c) is primarily due to overproduction of summer and fall precipitation. The summer season is also responsible for the overproduction of precipitation through the Bay of Bengal and Gulf of Thailand, suggesting an overactive summer monsoon in the MPAS-A simulations. Another notable difference between MPAS-A and TRMM annual average precipitation is the westward shift of the heaviest precipitation along the ITCZ in the Atlantic basin. This shift in precipitation is likely related to a westward shift of the summertime African easterly jet (AEJ; not shown).

Figure 6 Average total annual precipitation (mm d⁻¹) for (a) present-day MPAS-A simulations, (b) TRMM 19-year climatology, and (c) the model bias (MPAS-A minus TRMM). MPAS-A output was linearly interpolated to the TRMM grid for point-to-point comparison. Pattern correlation coefficient is reported in the top right of panel (b). Stippling in panel (c) indicates locations where the MPAS-A 10-year mean exceeds the range computed from 100 random samples of 10-year means from TRMM by more than 5 % (∼1.2 mm d⁻¹).

4 Model climate: tropical cyclones

Tropical cyclones epitomize the high-impact weather phenomena that our simulations are designed to address; therefore, we consider tropical cyclones as an appropriate exemplar of high-impact weather systems, in order to explore the usefulness of our simulations for examining the effects of climate change on such phenomena. To analyze how TCs appear in our model, we track simulated Northern Hemisphere tropical cyclones using the TempestExtremes objective, feature-based tracking algorithm (Ullrich and Zarzycki, 2017; Zarzycki and Ullrich, 2017). TCs are initially detected as minima in SLP, and then retained as candidate cyclone centers if certain criteria are met. Here, we require that TCs must have a 2 hPa closed SLP contour within 2^∘ of the storm center and a 300–500 hPa geopotential thickness maximum within 6^∘ of the storm center to ensure the presence of a warm core. Additionally, TCs must not travel more than 6^∘ within a 6 h period, must have a lifetime of at least 2 d, must be located over water for at least 12 h, must have at least 2 d of 10 m wind speed of at least 14 m s⁻¹, and are required to have a genesis latitude south of 45^∘ N. Trajectories that end and begin within 12 h of each other are merged together to prevent broken tracks from being counted twice. Once TC tracks have been obtained, TCs are separated into basins (Fig. 1) based on their genesis location.

4.1 Strength

We compare simulated TC characteristics to the International Best Track Archive for Climate Stewardship (IBTrACS; Knapp et al., 2010). Only the IBTrACS values for the 10 simulated years (Table 1) are considered for comparison. Consistent with similar studies (e.g., Murakami et al., 2015; Roberts et al., 2015; Yamada et al., 2017), model storms are generally weaker than observed in terms of maximum 10 m wind speed; several simulated storms do, however, attain a minimum SLP of less than 900 hPa (Fig. 7). Therefore, as in Roberts et al. (2015), storm intensity for the simulated TCs is measured by the minimum lifetime SLP of the storm in addition to maximum 10 m wind speed as defined by the Saffir–Simpson scale. Using the minimum SLP, categories are defined as >994 hPa for tropical storms (TS_p) and 980–994, 965–979, 945–964, 920–944, and <920 hPa for category (Cat_p) 1–5 equivalent tropical cyclones, respectively (Roberts et al., 2015). The subscript p is used to discriminate the SLP-based categories from those defined by the Saffir–Simpson wind speed thresholds. For IBTrACS, the maximum 10 m wind speed and minimum SLP across all reporting centers are used as the observed storm intensity for categorization.

Figure 7 Scatter plot of maximum 10 m wind speed (kts) versus minimum SLP (hPa) for IBTrACS (black) and present-day MPAS-A (grey) Northern Hemisphere TCs. The lines of best fit for each (IBTrACS in black and MPAS-A in red) were computed using a second-order polynomial. The wind (SLP) category thresholds are indicated by the vertical (horizontal) lines.

Figure 8 Average number of TCs over the 10 simulation years for the Northern Hemisphere and each Northern Hemisphere basin. Columns are colored by intensity categories based on (a) minimum lifetime SLP and (b) maximum lifetime 10 m wind speed. The bottom color represents intensities of tropical storm strength or less for IBTrACS and MPAS-A in the first and second columns, respectively. Categories 1–5 are shaded for both datasets according to the legend. The error bars indicate the interannual standard deviation. The number of TCs for IBTrACS varies based on strength metric due to the lack of SLP records for a select number of storms.

Figure 8 shows the average TC frequency over the 10 simulation years for the Northern Hemisphere as a whole, in addition to each basin. Our MPAS-A simulations generate excess TC activity in the Northern Hemisphere, primarily due to the overactivity in the western North Pacific basin. Simulated TC frequencies for the North Atlantic, eastern North Pacific, and Northern Indian basins are within the observed range. Across all basins, when categorizing TCs by minimum SLP (Fig. 8a), MPAS-A generally underestimates the number of weak systems (those with strengths less than Cat_p1), and overestimates the number of Cat_p1 and Cat_p2 storms. The frequency Cat_p3 TCs and stronger, on the other hand, are simulated reasonably well. With regard to TC categorization by maximum 10 m wind speed (Fig. 8b), MPAS-A simulates the frequency of TS strength TCs quite well in all basins. Strong TCs (Cat4 and Cat5), however, are universally underestimated by the model in favor of Cat1–Cat3 TCs. In the future simulations, there is an increase in TC activity in both the North Atlantic and western North Pacific basins (not shown); further investigation into these future changes in TCs will be the subject of a future paper.

4.2 Location

Spatially, the model-simulated TC track density compares reasonably well with observations; the pattern correlation coefficient is about 0.7 (Fig. 9). The most prominent difference is the lack of TC activity in the eastern portion of the North Atlantic basin, which is common among several similar modeling studies (e.g., Bell et al., 2013; Strachan et al., 2013; Small et al., 2014; Roberts et al., 2015). TC genesis in this region typically occurs during August and September (Kossin et al., 2010; Daloz et al., 2015); comparison between the simulated atmosphere and ERA-I monthly-averaged 850–200 hPa vertical wind shear for the 10 simulation years during these months shows a strong positive bias in the model over the North Atlantic development region that is likely a primary factor in this lack of TC generation (not shown). Additionally, the westward shift in the simulated summertime AEJ (not shown) is consistent with a westward shift in the wave accumulation zone, which is likely impacting the location of TC genesis in this area (Done et al., 2011). Figure 9 does not show a strong track density bias in the Gulf of Mexico. Roberts et al. (2015) note that a steady supply of vorticity in the Caribbean contributed to their overestimation of track density in this area; thus, it is possible that tracking TCs as SLP minima, rather than maxima in 850 hPa relative vorticity, helps alleviate this bias. Unlike previous studies (e.g., Small et al., 2014; Murakami et al., 2015; Roberts et al., 2015; Yamada et al., 2017), we do not find a positive track density bias in the central North Pacific; instead, we see a slight underrepresentation of TC activity in that area around ∼150^∘ W.

Figure 9 Track density (number of cyclone tracks per $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ area) over the 10 simulated years for (a) present-day MPAS-A simulations and (b) IBTrACS. Contours are shaded every one count.

4.3 Seasonal cycle

Aside from the underestimation in August and September (likely attributed to the lack of TC genesis in the eastern portion of the basin), MPAS-A simulates the present-day seasonal TC cycle for the North Atlantic reasonably well; TC activity increases during the spring and summer seasons and reaches a maximum in the fall (Fig. 10a). For the eastern North Pacific, MPAS-A produces too many storms in the springtime (April and May) and too few storms during the summer months (Fig. 10b). As defined by Camargo et al. (2008), cluster-2-type eastern Pacific TCs form off the coast of Mexico, travel towards the northwest along the coastline, and have a bimodal seasonal distribution with peaks in late spring/early summer and early fall, similar to the modeled cycle in Fig. 10b. Compared to the ERA-I 10-year climatology, enhanced westerlies at the 500 and 850 hPa levels between 0–20^∘ N in the eastern portion of the eastern Pacific basin for April and May (not shown) suggest that our simulations may be in a regime more conducive to these cluster 2 storms. For the western North Pacific (Fig. 10c), MPAS-A correctly simulates the fall peak in TC activity; there is, however, a secondary peak in April that does not match observations. Although there is a general overestimation of storm activity in the northern Indian basin, the model does replicate the shape of the seasonal cycle with both the early summer and mid-fall peaks represented, albeit the fall peak occurs 1 month earlier than observed (Fig. 10d).

Figure 10 Monthly-average TC frequency over the 10 simulated years for the (a) North Atlantic, (b) eastern North Pacific, (c) western North Pacific, and (d) northern Indian basins. The frequencies for IBTrACS (MPAS-A simulations) are shown in the grey (red) bars.

Figure 11 (a, b) Average annual 2 m temperature difference (K; future minus current) for the (a) MPAS-A simulations and (b) CMIP5 GCM ensemble mean. (c, d) Average annual difference cross-section of zonal mean temperature (K; future minus current) for the (c) MPAS-A simulations and (d) CMIP5 GCM ensemble mean. Contours are shaded every 1 K in all panels.

5 Climate change representation

To ensure that our simulations are useful in studying climate change effects on weather phenomena, we compare temperature change fields with large-scale warming patterns generated by a subset of IPCC GCMs. Previous theoretical and modeling studies demonstrate that the Arctic region will continue to warm at a faster rate than the rest of the globe in response to an increase in greenhouse gases (IPCC, 2013, Sect. 12.4.2.2). This polar amplification effect is captured in our simulations with portions of the Arctic experiencing temperature changes greater than 16 K compared to differences ≤10 K elsewhere (Fig. 11a). However, we note that Arctic temperatures in our present-day simulations were colder compared to ERA-I, while the future simulations resulted in Arctic temperatures comparable to those produced by GCMs (not shown). As a result, the MPAS-A simulations produce a larger magnitude of warming in the Arctic compared to the GCM ensemble (Fig. 11a–b). Another result consistent with theory and previous modeling studies is the development of a warming maximum in the tropical upper troposphere (IPCC, 2013, Sect. 12.4.3.2). This area of warming, which occurs between the ∼400 and ∼150 hPa levels and which maximizes around the 250 hPa level (IPCC, 2013, Sect. 12.4.3.2), has been shown to partially mitigate projected increases in TC intensity associated with warming (e.g., Knutson and Tuleya, 1999; Shen et al., 2000; Hill and Lackmann, 2011). As shown in Fig. 11c, this warming signature is replicated in the MPAS-A simulations.

6 Summary and conclusions

We present a novel set of model simulations produced using MPAS-A in current and future environments that is designed to maximize our ability to analyze changes in high-impact weather systems; such changes will be reported in future studies. Our use of a global model eliminates the lateral boundary constraints of regional models, while inclusion of high-resolution, analyzed SSTs preserves realistic SST gradients throughout the duration of the simulations. Furthermore, a grid length of 15 km offers an advantage over coarser modeling studies to better represent the atmospheric mesoscale. The future climate simulations employ a technique that combines methods associated with PGW and time-slice experiments; this allows for the inclusion of high-resolution SSTs, plausible future sea ice fields, and seamless simulation of non-consecutive years without excessive spin-up time. While the primary purpose of our simulations is to study climate change effects on Northern Hemisphere high-impact weather, to achieve this it is necessary to first evaluate the model climate in regard to present-day large-scale circulations as well as the large-scale responses to warming in our climate change experiments; reasonable representation of these aspects is essential to justify moving forward to investigate smaller-scale, high-impact phenomena.

Key results from these simulations include the ability of MPAS-A to reproduce Northern Hemisphere wintertime midlatitude storm tracks (Fig. 4) along with semi-permanent maritime SLP and upper-tropospheric jet features (Fig. 5). Tropical characteristics, such as precipitation along the ITCZ in the equatorial Pacific, are also well simulated (r∼0.95), although the ITCZ representation in the Atlantic does not compare as favorably to observations (Fig. 6). In regard to TC strength, MPAS-A is able to produce several tropical cyclones of Cat4 strength, as defined by traditional maximum 10 m wind speed thresholds of the Saffir–Simpson scale (Figs. 7 and 8b). Categorizing TCs using the minimum SLP thresholds of Roberts et al. (2015), on the other hand, shows simulated TCs across the full intensity spectrum, including Cat_p5 storms (Figs. 7 and 8a).

While MPAS-A overestimates TC activity in the western North Pacific, TC frequency in other Northern Hemisphere basins is within the range of observations (Fig. 8). The largest discrepancy in the simulated spatial distribution of TCs is the lack of TC genesis in the eastern North Atlantic (Fig. 9), likely due to a positive bias in vertical wind shear (not shown). Otherwise, TC density patterns match observations reasonably well (r∼0.7). Additionally, with the exception of the eastern North Pacific, the seasonal cycles for the Northern Hemisphere basins are well reproduced (Fig. 10). Last, our future simulations replicate two key warming signatures produced by GCMs: Arctic amplification and the warming maximum in the tropical upper troposphere (Fig. 11).

With our modeling approach, we strive to contribute to the intersection of weather and climate modeling, and aim to fill a gap between GCMs, which are unable to simulate small-scale weather phenomena, and high-resolution limited-area models, which are constrained by lateral boundaries, to provide the possibility of studying high-impact events in a consistent global context. We anticipate these simulations, in conjunction with similar efforts, will have great value in projecting and understanding changes in high-impact weather phenomena for which dynamics on sub-synoptic scales are important. Beyond the tropical cyclones described here, this could include flooding rains and damaging winds associated with extratropical cyclones, flooding monsoon rains, and localized droughts and heat waves. Research involving these simulations is currently underway, investigating climate change effects on the following phenomena:

• extratropical transition of TCs,

• TC seasonality,

• midlatitude precipitation extremes and windstorms embedded in extratropical cyclones, and

• persistent anomalies and blocking.

Many more aspects of these simulations, however, remain to be explored. Therefore, we make the simulation output available to the research community as detailed in the code and data availability section in the hope that it will be useful to the broader scientific community for studying various meteorological phenomena, as well as for conducting model comparison studies.