SPEAR – the next generation GFDL modeling system for seasonal to multidecadal prediction and projection

©2020 American Geophysical Union. All rights reserved. Delworth Thomas, L. (Orcid ID: 0000-0003-4865-5391) Cooke William, F. (Orcid ID: 0000-0003-4550-3210) Adcroft Alistair (Orcid ID: 0000-0001-9413-1017) Bushuk Mitchell (Orcid ID: 0000-0002-0063-1465) Chen Jan-Huey (Orcid ID: 0000-0002-1181-8110) Ginoux Paul (Orcid ID: 0000-0003-3642-2988) Harris Lucas (Orcid ID: 0000-0001-6072-8624) Harrison Matt (Orcid ID: 0000-0001-7991-454X) Johnson Nathaniel (Orcid ID: 0000-0003-4906-178X) Kapnick Sarah, B (Orcid ID: 0000-0003-0979-3070) Malyshev Sergey (Orcid ID: 0000-0001-6259-1043) Naik Vaishali (Orcid ID: 0000-0002-2254-1700) Paynter David (Orcid ID: 0000-0002-7092-241X) Schwarzkopf M. Daniel, Daniel (Orcid ID: 0000-0003-1381-1219) Underwood Seth (Orcid ID: 0000-0002-8395-9609) Wittenberg Andrew, T. (Orcid ID: 0000-0003-1680-8963) Xiang Baoqiang (Orcid ID: 0000-0001-5774-7437) Yang Xiaosong (Orcid ID: 0000-0001-5601-9710) Zhang Liping (Orcid ID: 0000-0003-3548-321X) Zhao Ming (Orcid ID: 0000-0003-4996-7821)

with available computational resources. For computational speed SPEAR uses a coarse ocean resolution of approximately 1.0 o (with tropical refinement). SPEAR can use differing atmospheric horizontal resolutions ranging from 1 o to 0.25 o . The higher atmospheric resolution facilitates improved simulation of regional climate and extremes. SPEAR is built from the same components as the GFDL CM4 and ESM4 models, but with design choices geared toward seasonal to multidecadal physical climate prediction and projection. We document simulation characteristics for the time-mean climate, aspects of internal variability, and the response to both idealized and realistic radiative forcing change. We describe in greater detail one focus of the model development process that was motivated by the importance of the Southern Ocean to the global climate system. We present sensitivity tests that document the influence of the Antarctic surface heat budget on Southern Ocean ventilation and deep global ocean circulation. These findings were also useful in the development processes for the GFDL CM4 and ESM4 models. ©2020 American Geophysical Union. All rights reserved.

Plain Language Summary
In this paper we describe the development and simulation characteristics of a new climate model that will be used for seasonal to multidecadal climate prediction and projection. The model combines a set of newly developed components that simulate the ocean, atmosphere, land and sea ice. We document this model by assessing its performance in simulating the current climate, and by showing the model's response to changing greenhouse gases and aerosols over the 20 th and 21 st centuries. We also show results from a set of sensitivity experiments that were an important part of the model development process. These sensitivity tests explore connections between the surface energy balance over Antarctica and the circulation of the deep ocean.
In this manuscript we describe the development process and simulation characteristics of a new modeling system that will form the backbone for the next generation real-time prediction and projection efforts at the Geophysical Fluid Dynamics Laboratory (GFDL), extending seamlessly from seasonal to multidecadal time scales. Real-time experimental seasonal predictions have been conducted since the early 1990s at GFDL, with real-time predictions routinely distributed to institutions such as the International Research Institute for Climate and Society at Columbia University. Beginning in 2015 this ongoing activity contributed model predictions to the newly formed North American Multimodel Ensemble (Kirtman et al., 2014) seasonal prediction activity. Similarly, real-time decadal predictions have been performed routinely at GFDL since 2011 (Yang et al., 2013), with model predictions currently being made available through the UK Met Office as a WMO Lead Centre for Interannual to Decadal Prediction (https://www.metoffice.gov.uk/research/climate/seasonal-to-decadal/long-range/wmolcadcp). These activities are part of a seamless seasonal to centennial climate prediction and projection effort at GFDL, as the same models that are used for seasonal and decadal prediction have also been used for century-scale climate change projections (Held et al., 2010(Held et al., ,2019Knutson et al., 2006;Merlis et al., 2014).
The two coupled climate models currently used for the real-time experimental predictions are CM2.1  and FLOR (Vecchi et al., 2014). FLOR is derived from CM2.1 but uses a higher resolution atmospheric component (50 km horizontal grid spacing in FLOR versus 200 km in CM2.1). Both coupled systems use an ocean model with a horizontal resolution of approximately 1 o , with refined meridional resolution in the Tropics. The seasonal to decadal predictions are initialized using the output of an Ensemble ©2020 American Geophysical Union. All rights reserved.
Coupled Data Assimilation (ECDA) system (S. Zhang et al., 2007). The seasonal predictions are computed in real-time each month, using all data from the previous month.
The fundamentals of the CM2.1 and FLOR models were developed more than 10 years ago. Since then there has been substantial component model development and improvement at GFDL. These new developments provide the rationale and basis for creating a new system for seasonal to multidecadal predictions and projections at GFDL that takes advantage of these advances. These developments include the adoption of a new dynamical core in the atmosphere with regional nesting capabilities (Harris & Lin, 2013), revised atmospheric physics (Zhao et al., 2018a(Zhao et al., , 2018b, the new MOM6 ocean code and SIS2 sea ice code  , and a substantially enhanced land model. These new developments will enable new capabilities and advances. For example, the new dynamical core will allow for the use of variable resolution grids, thereby allowing future versions of SPEAR to use high spatial resolution in focused regions for the simulation and prediction of critical phenomena, such as hurricanes or other extremes. The MOM6 ocean code includes many improvements, such as the use of hybrid vertical coordinates, that may assist in the simulation of important oceanic processes potentially relevant to decadal prediction, such as overflows. A much improved representation of land hydrology will enable the simulation and prediction of extremes in river flows. This manuscript describes the development process and simulation characteristics of this new seasonal to multidecadal prediction and projection system, which we call SPEAR (Seamless System for Prediction and EArth System Research). The SPEAR models share many components with the recently developed GFDL CM4 (Held et al., 2019) model, but with configuration and physical parameterization choices in SPEAR geared toward physical climate prediction and projection on seasonal to multidecadal time scales. Table 1 places SPEAR within the context of a series of models developed at GFDL over the last 15 years.
In section 2 we describe the model components and the simulations conducted. We present in section 3 (and accompanying online supporting information) a suite of diagnostics to illustrate model characteristics in control simulations, as well as the model response to idealized and realistic radiative perturbations. One aspect of the model development process that deserves special attention, which we focus on that in section 4, has to do with the role of the Southern Ocean (SO) in the global climate system, and its influence on model bias and drift. We illustrate aspects of the development process that focus on simulating the Southern Ocean and its global-scale climatic impacts through the design and analysis of sensitivity experiments that illustrate connections between continental climate over Antarctica, the Southern Ocean, and model bias and drift. A summary is presented in section 5.

Model Description
We use a newly coupled ocean-atmosphere-land-sea ice modeling system (SPEAR) that was developed as a seamless climate prediction system capable of making predictions and projections from seasonal to centennial time scales. The Earth's climate system is seamless and unified, with short-term predictable variations, such as El Niño-Southern Oscillation (ENSO) or the Madden-Julian Oscillation (MJO), occurring at the same time as the system responds to changing greenhouse gases, aerosols, and other forms of radiative forcing changes. Further, changing radiative forcing can in turn alter the characteristics of internal variability of the climate system (Fang et al., 2014;Zhang & Delworth, 2016b).
Modeling systems need to be designed to predict the behavior of the climate system as an initial value problem, the response of the climate system to radiative forcing changes, and the influence of radiative forcing changes on the characteristics of internal variability. This latest model is part of a tightly coupled new suite of models developed recently at GFDL (see Table   1). These models share the vast majority of their code and components but are designed with specific foci in mind. The SPEAR modeling system was designed with seasonal to multidecadal prediction and projection as its focus, with an emphasis on computational efficiency to facilitate large ensembles needed for predictions and projections.

Atmosphere/land formulation
The physics in the atmospheric and land components of SPEAR are identical to those of the recently developed GFDL AM4-LM4 model (Zhao et al., 2018a(Zhao et al., , 2018b, with one exception: the near-infrared isotropic reflectance parameter for cold snow over glacial surfaces has been increased from 0.58 to 0.68 (see Milly et al (2014) for details on the formulation for snowpack in the land component). The increased albedo, as discussed in more detail in section 4, both improves the realism of the simulated albedo over Antarctica and improves aspects of the circulation of the World Ocean related to deepwater formation in the SO and global ocean temperature drift.
We use two versions of SPEAR that differ in the horizontal resolution for the atmosphere and land. SPEAR_LO uses an atmosphere/land resolution of approximately 100 km, and SPEAR_MED uses an atmosphere/land resolution of approximately 50 km. The 100km resolution of SPEAR_LO is identical to that in Zhao et al (2018), with the exception of the albedo change noted above. The physics in the 50-km version of SPEAR are identical to those in the 100-km version. The dynamical core differs only in the timesteps required for the difference in resolution. The atmospheric component uses 33 levels in the vertical. We present results in Sections 2 and 3 (and online supplemental information) using both resolutions. For a more detailed discussion of the physics of the atmosphere and land components, as well as model performance when forced by observed sea surface temperatures (SST), the reader is referred to the detailed results in Zhao et al. (2018). We present results below (and in the online Supporting Information) for model performance when coupled to a dynamic ocean. Both SPEAR_LO and SPEAR_MED are coupled to the same 1 o ocean model as described in the next section.

Ocean/sea ice formulation
The ocean and sea ice components are constructed from the new MOM6 ocean code (https://github.com/NOAA-GFDL/MOM6) and SIS2 sea ice code (https://github.com/NOAA-GFDL/SIS2). The ocean and sea ice components use a nominal horizontal resolution of 1 o with refinement to 1/3 o in the meridional direction in the Tropics.
Since the model grid maintains an aspect ratio close to 1, the effective resolution increases poleward with the convergence of the meridians. The zonal grid spacing reduces from 110 km in the deep tropics to 46 km at 65 o latitude. The meridional grid spacing is 40 km in the deep tropics (this equatorial refinement is important for representing tropical variability), increases to 105 km in the subtropics, and decreases to 46 km at 65 o latitude in each hemisphere. In order to maintain approximately uniform grid spacing over the Arctic, a tripolar grid is used with poles over northern Eurasia, northern Canada, and Antarctica (Murray, 1996). North of 65 o N there is a non-uniform grid with cells of approximately equal area (approximately 46 km by 46 km) covering the polar cap. The ocean model has 75 layers in the vertical, with layer thickness as fine as 2 m near the surface, including 30 layers in the top 100 m. The vertical grid spacing becomes larger with depth, reaching 250 m below 5000m.
The ocean resolution employed in SPEAR is relatively coarse. This is a deliberate choice to facilitate the rapid and efficient execution of ensembles of simulations needed for climate prediction given available computer resources. Based on past experiences with FLOR we also decided to employ higher resolution in the atmosphere (as fine as 50 km with planned refinement to 25 km) in order to simulate important phenomena such as tropical storms, other extremes, and regional hydroclmate. We plan to explore future versions of SPEAR with enhanced ocean resolution as resources permit.
The physics used in the ocean and sea ice components is very similar to that documented in Adcroft et al., 2019 (hereafter called OM4, the ocean component of CM4), including the use of a hybrid vertical coordinate in the ocean. The hybrid vertical coordinate is a function of height in the upper ocean, transitioning to isopycnal layers in the interior (see Adcroft et al, 2019, for further details). The use of the hybrid vertical coordinate provides a more realistic solution than model versions in our development process that used a purely level based coordinate system; improvements include a reduced surface temperature bias and less subsurface temperature drift. Both models (SPEAR and OM4/CM4) use the same 75 vertical levels. There are, however, a few key differences reflecting the differences in horizontal resolution between OM4, which uses a 0.25 o horizontal grid spacing, versus the 1 o grid spacing used in the SPEAR model. First, the SPEAR model is not able to resolve oceanic mesoscale eddies, and thus uses a parameterization of their effects. Spatially varying model lateral diffusivities are internally diagnosed based on a mesoscale energy budget. The production and dissipation of sub-grid scale mesoscale kinetic energy balanced by the release of resolved available potential energy is parameterized by an eddy overturning streamfunction (Marshall & Adcroft, 2010;Jansen et al., 2015). Additionally, predicted eddy velocity scales are combined with an eddy mixing length scale to calculate tracer mixing rates along neutral density surfaces. Second, similar to OM4, the SPEAR ocean component uses a parameterization of the effects of submesoscale eddies  . Relative to OM4, different parameter coefficients were used in SPEAR based upon the results of sensitivity tests and the impact on SST bias. Third, the MOM6 code has the capability to represent flow through channels that are narrower than the grid spacing of the model. The locations where this functionality is needed are dependent on the model grid, and this was set up for SPEAR based on the grid resolution at 1 o . These include regions with narrow channels such as occur in the Indonesian Throughflow, the Gibraltar Straits, and the Bosporus strait.
This parameterization was tuned to produce flow through the channels similar in total transport to observed. Fourth, we found empirically an improvement in some aspects of the model solution by adding extra horizontal viscosity poleward of 50 o in each hemisphere.
Highly varying topography and a declining Rossby radius of deformation at high latitudes are challenges for a coarser resolution model, and our experience was that the extra viscosity provided a better solution.

Observational data
We use observational data sets for model assessment. For sea surface temperature (SST) we use NOAA Extended Reconstructed Sea Surface Temperature version 4 (ERSSTv4) data (Huang et al., 2015) provided by NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, from their Web site at https://www.esrl.noaa.gov/psd/. For surface air (surface air temperature over land and sea surface temperature over the ocean) we used two data sets. The first (CRUTEM4) was developed at the Climatic Research Unit (Jones et al., 2012) and downloaded April 26, 2019, from NOAA's Earth System Research Laboratory (https://www.esrl.noaa.gov/psd/data/gridded/data.crutem4.html). The second (GISTEMP v4) is from NASA's Goddard Institute for Space Sciences (GISS) (Lennsen et al., 2019), and was downloaded on July 30, 2019 from https://data.giss.nasa.gov/gistemp/. For sea ice extent we used observational data based on passive microwave satellite information (Cavalieri et al., 1996) obtained from the National Snow and Ice Data Center derived from brightness temperature measured by satellite using the NASA Team algorithm. Atmospheric circulation ©2020 American Geophysical Union. All rights reserved.

Control simulations
We conduct extended control simulations with atmospheric composition fixed at levels representative of either calendar year 1850 (a "preindustrial" control) or calendar year 2010 (a "modern" control). In each case the model ocean is initialized from an observed late 20 th century ocean climatology (Boyer et al., 2013). We conduct control simulations using both SPEAR_LO and SPEAR_MED. As shown below, in the 1850 control simulations the atmosphere and ocean cool over the first century as the system adjusts to 1850 radiative forcing. The initial conditions come from a late 20 th century observed state with a positive net radiative imbalance (Allan et al., 2014), whereas 1850 radiative conditions have a radiative imbalance closer to zero. This difference in radiative forcing would be expected to lead to cooler conditions for the 1850 control simulations than observed in recent decades. We show key aspects of the model solutions to provide an overall sense of model fidelity in simulating  As with AM4 and CM4, we specify time series of stratospheric aerosol optical properties including contributions from both volcanoes and other natural and anthropogenic contributions (Thomason et al., 2018). For the period 2025-2100 we specify stratospheric aerosol properties as the time-mean used over the 1850-2014 period. During the period 2015-2024 we linearly interpolate from the historical values in 2014 to the time-mean values specified over the period 2025-2100. The contribution from continuously degassing and explosive volcanoes to tropospheric aerosols is treated in the same way as in AM3  and AM4 (Zhao et al., 2018a).
Tropospheric aerosols in SPEAR are simulated from emissions with a simplified chemical mechanism as implemented in the AM4 model (Zhao et al., 2018a). A brief description of the representation of aerosols is provided here with more details given by Zhao et al. (2018a). Mass distribution of five aerosol types including sulfate, dust, black carbon, organic carbon, and sea salt are simulated in SPEAR. The size distribution for all aerosols is prescribed as lognormal except for dust and sea salt, which are discretized into five size bins from 0.1 to 10 micro m radius. Their concentrations are calculated based on their emissions (and precursor emissions), chemical production for sulfate and secondary organics, dry and wet (rainout and washout) deposition, transport by advection, and dry and wet convection.
Chemical production of sulfate is from the oxidation of sulfur dioxide (SO2) in the gas phase by hydroxyl radical (OH), and in the aqueous phase by ozone and hydrogen peroxide (H2O2).
Since ozone chemistry is not explicitly represented in the model, the monthly climatology of ozone and other oxidants (OH, HO2, and NO3) used by sulfate chemistry is prescribed based on a 20-year simulation of AM3 (Naik et al., 2013) that applied a comprehensive mechanism of ozone chemistry (Horowitz et al., 2003). Emissions of natural aerosols including dust, sea- positive biases off the coast of Antarctica, also associated with difficulties in simulating cloud cover, and (c) negative biases over some tropical continental locations. The global mean, annual mean net radiation at the top of the atmosphere (averaged over years 1-100) is 0.80 Wm -2 in SPEAR_LO and 0.70 Wm -2 in SPEAR_MED, in good agreement with recent observational estimates (Johnson et al., 2016;Loeb et al., 2009Loeb et al., , 2012 . These positive imbalances reflect the effect of increased CO2 and other radiative forcing agents in calendar year 2010 relative to a longer-term steady state (the radiative balance associated with preindustrial conditions is discussed below). the Gulf Stream and its separation from the North American coast. This is a typical bias in models that do not explicitly resolve oceanic mesoscale eddies, see Figure 9.2 of Flato et al. (2013) . Positive biases off the west coast of Africa are associated with difficulties in simulating marine stratocumulus decks and regional atmospheric circulation, as reflected in the errors in the net radiation balance (see Figure 1). Both problems are long-standing biases in many climate models (see, for example, Figure  The annual mean precipitation is shown in Figure 3 for both observations and models. While the overall patterns are realistic, there is still a tendency for a double InterTropical Convergence Zone (ITCZ) in the tropical Pacific, along with excessive rainfall in the tropical Atlantic and the western tropical Indian Ocean. We examine in more detail the ability of the model to simulate regional precipitation patterns over North America in Figure 4. The patterns are realistic, but there is an overall wet bias, with more simulated precipitation than observed over North America. As expected, the finer resolution of SPEAR_MED provides more details in the simulation of orographically influenced precipitation over western North America, including regional scale dry zones in the lee of mountainous areas.
The simulation of zonal winds is very important in driving ocean circulations, and simulated winds are shown in Figure 5. The top panel shows the zonal mean of the nearsurface zonal wind stress over the Pacific Ocean sector; common to many models the zonal winds in the mid and higher latitudes of the Southern Hemisphere (SH) are shifted equatorward in the models relative to their locations in observations-based reanalyses (Dee et al., 2011;Rienecker et al., 2011), and see also Figure 9.19 of Flato et al. (2013). This bias is not as readily apparent in the Atlantic sector (middle panel), although simulated zonal winds in the middle and higher latitudes of the Northern Hemisphere (NH) also have an equatorward shift relative to observations. Over continental regions (bottom panel) the locations and amplitudes of the simulated winds are in good agreement with observational estimates. A further perspective on simulated atmospheric circulation is shown in Figure 6 which shows sea level pressure (SLP) for the months of Dec-Feb (DJF) over the Northern Hemisphere. Both models have a tendency for southward displacements of the Icelandic and Aleutian Lows relative to observations, consistent with the zonal wind stress shown previously. ©2020 American Geophysical Union. All rights reserved.
The seasonal cycles of sea ice extent in observations and the model simulations are shown in Figure 7. Simulations in the Northern Hemisphere are fairly realistic. There is somewhat too much winter sea ice extent, especially in the Greenland Sea and Barents Sea regions (not shown). In the Southern Hemisphere the model simulates a larger than observed seasonal cycle of sea ice, with excessive extent in winter and summer values below observed.
This bias may be related to a tendency for excessive absorbed shortwave radiation in summer, thereby reducing sea ice, or to a lack of dynamical thickening (ridging) processes in this model  .

Time-varying climate from internal variability
We also conduct 2000-year control simulations with atmospheric composition set at 1850 conditions. The lower concentrations of atmospheric greenhouse gases in 1850 relative to 2010 produces a simulated climate in 1850 that is both cooler, and closer to radiative balance, than the simulated 2010 climate. The long-term mean in global mean net radiation at the top of the atmosphere over years 501-1000 is 0.18 Wm -2 in SPEAR_LO, and 0.12 Wm -2 in SPEAR_MED. The smaller imbalance in SPEAR_MED appears to be related to differences in the cloud simulations between the two models and is a focus of current diagnostic analyses.
The models do undergo a long adjustment process, as the initial conditions from the ocean are representative of modern conditions, and the simulated climate system must move towards an equilibrium with atmospheric composition (and hence radiative forcing) consistent with the preindustrial period. We show in Figure 8a the time series of the net radiative balance at the top of the atmosphere for the 1850 control simulations of SPEAR_LO and SPEAR_MED. After an initial rapid adjustment, both models have a net positive bias that slowly declines over time, indicating that the climate system is gaining energy. This is a typical bias in many previous GFDL climate models. Even with 1850 radiative forcing (which should be less positive than today's climate) the ocean warms by moving heat from the near-surface layers into the ocean interior. Aspects of this process have been explored using an earlier coupled climate model with very high ocean resolution (0.1 o ) (Delworth et al., 2012) in which it was shown that an improved representation of the effects of ocean mesoscale eddies can ameliorate this problem. This persistent bias points to the importance of either explicitly resolving fine scale processes in the ocean, including mesoscale eddies, or developing improved parameterizations of their impacts for use in coarser resolution ocean models.
We show in Figure  We next examine the model's ability to simulate a variety of phenomena relevant for seasonal to decadal predictions, starting with the Madden-Julian Oscillation (Madden & Julian, 1971). As one measure of the ability of SPEAR to simulate the MJO, we show in On seasonal time scales the El Nino/Southern Oscillation (ENSO) is a dominant phenomenon. We show in Figure 10 the standard deviation of monthly sea surface temperature (SST) anomalies in the tropical Pacific. The models broadly capture the spatial structure of SST variability in the tropical Pacific, although both models appear to have variability that is somewhat larger than observed and that extends further westward in the models than in observations. The higher resolution SPEAR_MED model appears to have a slightly better simulation then SPEAR_LO by this metric. The time scale of ENSO variability is shown in Figure 11 by the spectra of SST averaged over the Niño 3 region the Pacific provide realistic simulations of the teleconnections, although there is a tendency for the Pacific and North American maximum teleconnections to be shifted westward relative to observations. The same is true for the previous generation GFDL prediction model (FLOR, see Table 1) as shown in Figure 12b. 2016) in the prediction of seasonal tropical storm and hurricane activity. We show in Figure   13 the ability of SPEAR_LO and SPEAR_MED to simulate the frequency of tropical storms and hurricanes in various basins around the globe compared to observations (Knapp et al., 2010) . While the spatial resolution of both models is insufficient to simulate major hurricanes, both models are able to simulate the broad scale distributions of tropical storms in most ocean basins, with SPEAR_MED able to simulate more intense storms than SPEAR_LO. The spatial distribution of simulated storms is also more realistic in SPEAR_MED than SPEAR_LO.
On decadal time scales the Pacific Decadal Oscillation (PDO, also referred to as Pacific Decadal Variability, PDV) is the dominant pattern of variability in the Pacific (Mantua et al., 1997). This low frequency variability is associated with a large-scale change in the climate system and marine ecosystems (Cayan et al., 2001;Deser et al., 2004). We show in Figure 14 the spatial pattern and power spectrum of the PDO in observations and SPEAR, along with the older generation FLOR model (see Table 1). The SST anomaly associated with the PDO in observations is characterized by a horseshoe-like pattern over the North Pacific, with SST anomalies of one sign in the western and central North Pacific and SST anomalies of the opposite sign to the north, south and east ( Figure 14a). The observed pattern shows its maximum over the central North Pacific along with a secondary maximum in the Kuroshio-Oyashio-Extension (KOE) region. The SPEAR models generally capture the observed horseshoe-like shape of the PDO over the North Pacific (Figure 14 a, c and e).
However, the SST anomalies in the KOE region are higher than that in the central North Pacific. This issue is very common in fully coupled climate models, which is also seen in our previous FLOR model (Figure 14g). The North Pacific horseshoe-like pattern associated with the PDO also shows a connection with the tropical Pacific Ocean (Figures 14a, c, e, g). The SST anomalies extend from the eastern North Pacific to the central equator, suggesting strong teleconnections between the extratropic and tropic oceans. Compared to observation, this connection is overall weaker in models. However, the SPEAR models have a better performance compared to the FLOR model and is closer to observation. This improvement is related to better simulations of wind, SST and low clouds mean states over the eastern North Pacific in the SPEAR models. year, is not captured by the GFDL models. We also examined the impact of the PDO on North American precipitation. Consistent with observation and the FLOR model (Zhang & Delworth 2015), the PDO warm phase corresponds to wet conditions in the southwestern United States in both warm and cold seasons and vice versa. This is because the anomalous southerly associated with the Pacific-North-American (PNA) teleconnection during the PDO warm phase enhances moisture transport from south to north and therefore leads to an excess of precipitation and vice versa.
One of the prime processes leading to decadal predictability is the Atlantic Meridional Overturning Circulation (AMOC) and its slow decadal variations. We show in Figure 15 the AMOC for both SPEAR_LO and SPEAR_MED, calculated in two ways. In the first method the volume flow on height surfaces in the ocean (depth space) is used for computing the overturning streamfunction. In the second method the volume flow in density layers is used to compute the streamfunction. The AMOC calculated in depth space can be problematic, as northward and southward flowing water parcels at the same depth would cancel in the calculation of a net meridional flow. In reality those water parcels could have very different temperature and density characteristics so that there is a net poleward heat transport. The differences between the AMOC calculated in depth space and density space are most apparent at higher latitudes of the North Atlantic, where southward flowing cold denser water in the western Atlantic occurs at the same depth as northward flowing warm (light) waters in the eastern Atlantic, leading to a significant transport in density space but minimal transport calculated in depth space.

Ensembles of climate change simulations over the period 1851 to 2100
Although not the primary focus of this paper, we briefly document the simulation characteristics of the 30-member ensembles of the SPEAR_LO and SPEAR_MED models when driven by historical and projected changes in radiative forcing over the period 1851-2100 (for SPEAR_LO) and 1921-2100 (for SPEAR_MED). The design and radiative forcing for these simulations were discussed in section 2.4.3. These ensembles will be analyzed in greater detail in forthcoming papers, but we document here some basic aspects of their response to changing radiative forcing. Figure 17 are the time series of global mean surface air temperature from SPEAR_LO and SPEAR_MED, along with two observational estimates (CRUTEM4, 1850and GISTEMP v4, 1880. We show the 30-member ensemble mean from SPEAR_LO and SPEAR_MED as the thick solid red and blue lines respectively. The range of temperatures simulated across the ensemble members is denoted by the tan shading for SPEAR_LO, and the light blue shading for SPEAR_MED. The models reproduce many of the features of the observed record, including trends over the last 80 years. The models do not reproduce the early 20 th century warming seen in the observations. The global mean surface air temperature increase by the year 2100, relative to 1961-1990, is 4.54K for SPEAR_LO and 4.78K for SPEAR_MED when forced with the SSP5-85 scenario  of future radiative forcing.

Shown in
The patterns of changes in surface air temperature and annual mean precipitation are shown in Figure 18, calculated as the mean over the period 2041-2060 minus the period . Consistent with many previous model projections , warming is larger in the Northern Hemisphere than in the Southern Hemisphere (Fig 18a and 18b), with generally more warming over continental regions than oceanic regions. The muted warming over India and sub-Saharan Africa is related to increased rainfall (Figure 18c and 18d) and moister land surfaces, allowing any increase in the surface energy balance to be preferentially balanced by latent heat fluxes. The cooling south of Greenland is associated with a weakening of the Atlantic Meridional Overturning Circulation (AMOC, shown below in Figure 19) and associated reduction in oceanic meridional heat transport. The relative lack of warming in the Ross Sea and Southern Ocean is due to relatively strong vertical mixing in the ocean, leading to an effectively larger heat capacity and slower rate of temperature increase.
The precipitation changes (Figures 18c and 18d) have many aspects that have been seen in many previous models using earlier versions of future radiative forcing scenarios.  . The increasing aerosols preferentially cool the Northern Hemisphere and reduce poleward water vapor transport, thereby cooling the near-surface waters in the higher latitudes of the North Atlantic and increasing salinity due to reduced poleward water vapor transport and runoff into the ocean. The increased upper ocean density leads to enhanced deepwater formation and an increase in the AMOC. The long-term decline after the year 2000 is due to both declining atmospheric aerosols and increasing greenhouse gases; together these two factors warm and freshen the near-surface waters, thus increasing vertical stability and reducing deepwater formation and the AMOC.
The time series of the areal extent of Arctic sea ice in September is shown in Figure   20 for SPEAR_LO and SPEAR_MED, as well as the observations. The observed values lie within the model spread for SPEAR_LO, and partially within the spread for SPEAR_MED (the mean extent for SPEAR_MED is slightly larger than SPEAR_LO and the observations).
Trends of sea ice extent in recent decades are well simulated by both SPEAR_LO and SPEAR_MED. During the period 2007-2012 there appeared to be a more rapid decline of Arctic sea ice extent, but the observed trend over the last 5 years has moderated, so that the overall trend in both models agrees well with the observations. The first year in which the simulated Arctic in SPEAR_LO has ice-free conditions in September (defined as areal extent less than one million square kilometers) ranges from 2039 to 2060, with slightly later dates in SPEAR_MED due to its small positive bias in extent over the instrumental period.
The above analysis represents a brief assessment of these ensembles and the performance of the SPEAR models in simulating climate change from the 19 th through the 21 st centuries. The output from these large ensembles will be examined in much greater detail in future studies, and output from these simulations will be made publicly available. In particular, the 50 km resolution of SPEAR_MED makes this ensemble particularly useful for studies of regional hydroclimate and extremes.

Issues in simulating the Southern Ocean
As the emission of CO2 from industrial activity into the atmosphere continues, the amount of atmospheric CO2 that enters the SO and subsequently moves to depth is an important factor determining how much CO2 remains in the atmosphere, with important implications for the characteristics of future climate change (Frölicher et al., 2015;Landschützer et al., 2015;Shi et al., 2018). Due to relatively low vertical stratification, nearsurface water in the SO can move relatively easily to deeper layers of the ocean, thus facilitating the exchange of heat, fresh water, and CO2 between near-surface waters and the deep ocean. The characteristics of the upper ocean, including its variability, are tightly linked to the atmosphere via strong air-sea fluxes of heat, water and momentum. In addition, model simulations suggest that the SO is characterized by long time scales of variability that may be predictable (Latif et al., 2017;, and that the characteristics of the variability depend on the model base state. Further, characteristics of the SO are a crucial factor influencing the melting rate of the seaward edge of land-based ice sheets, with substantial implications for sea level change. Simulating a realistic SO is thus an important metric for climate model development. We characterize some aspects of circulation in the SO, including deepwater formation and ventilation, through the use of an ideal age tracer in the ocean. At the start of a simulation this tracer is set to zero everywhere in the world ocean. When a parcel of water is not in contact with the atmosphere the age tracer grows linearly with time (over 1 year the age tracer increases by 1 year); however, whenever a water parcel moves back to the surface and is in contact with the atmosphere the age tracer is once again set to zero. The age tracer is advected and mixed in the interior ocean similarly to a tracer such as salinity but has no influence on density or any other field and thus has no feedback on the ocean circulation.
Relatively low values of age indicate regions where water parcels have been in recent contact with the atmosphere. Vigorous AABW formation and ventilation should be characterized by low ages in the deeper layers of the SO, as water that has been recently at the surface moves to great depth. We show in Figure  Atlantic from the subpolar North Atlantic to the South Atlantic. In the SO we see extensive areas with age values less than 200 years, indicating substantial regions with recently ventilated water. As we move to the deeper ocean below 4000m (bottom row) we see the continued presence of low age values in the SO, indicative that recently ventilated water is sinking into the deepest layers of the SO. In contrast, in the North Atlantic the age values below 4000m are considerably larger/older than they were at 2000 m depth, indicating that the NADW formed in the subpolar gyre does not sink to the deepest layers of the North Atlantic. We also note that age in the deep North Pacific is smaller/younger than at 2000 m, indicating the propagation of younger water masses resulting from AABW formation into the deepest layers of the North Pacific.
The age values shown above were for the final version of SPEAR after a number of sensitivity tests related to SO ventilation. In early stages of the SPEAR model development process there was little deepwater formation and ventilation in the SO, leading to very low rates of AABW production, likely inconsistent with observations as inferred from oceanic tracers. The model has a warm fresh bias in the surface layers of the SO (associated with excess net shortwave radiation at the surface), thereby leading to vertical stratification that was larger than observed and inhibiting vertical exchange. Heat is transported in the subsurface ocean from the Atlantic Ocean into the SO; without sufficient vertical exchange this heat accumulates in the subsurface ocean. After sufficient time, this subsurface accumulation of heat in the model tended to destabilize the water column, eventually leading to abrupt releases of heat from the ocean through large, decadal-scale convective events that appeared inconsistent with the observational record.
In the model development process many sensitivity studies were done examining multiple processes related to SO deepwater formation, including tidal and other vertical mixing processes in the ocean. We found a crucial dependence of SO deepwater formation on the surface energy balance over Antarctica. The initial value specified for the near-infrared isotropic reflectance parameter for cold snow over glacial surfaces was 0.58, and resulted in albedos that were somewhat lower (~0.73) than suggested by observational estimates (Wendler & Kelley, 1988;Pirazzini, 2004;Laine, 2008;Wang & Zender, 2011) which range from 0.81 to 0.85. Increasing the reflectance parameter to 0.68 in the final version of SPEAR produced albedo values over Antarctica of ~0.81, in better agreement with the limited observational datasets. We found that these adjustments had a substantial impact on SO deepwater formation. We illustrate this aspect of the model development process by conducting a pair of idealized sensitivity studies in which we adjust this parameter over Antarctica and show how this impacted not only SO deepwater formation but characteristics of the global ocean and heat balance. Given the relative similarities between SPEAR_LO and SPEAR_MED we focus on results from SPEAR_LO.

Sensitivity tests
We conduct perturbation experiments using SPEAR_LO in which we modify the surface energy balance over Antarctica so that the air flowing off the continent is systematically cooled or warmed relative to a control simulation. In these experiments we alter the albedo of snow on glacial surfaces over Antarctica. For experiment "BRIGHT" we increase the near-infrared isotropic reflectance parameter (the main determinant of nearinfrared albedo) for cold snow over glacial surfaces south of 60 o S from 0.68 to 0.83. This is done only for continental locations. Similarly, for experiment "DARK", we reduce the value of this parameter from 0.68 to 0.53. In "CONTROL" we use a value of 0.68, as in SPEAR_LO. Both cases have a parameter perturbation of 0.15, but of differing sign. We anticipate that: (a) the BRIGHT (DARK) simulation will reflect more (less) radiation, and therefore will be colder (warmer) than the control simulation over Antarctica, and (b) the air flowing off the Antarctic continent in BRIGHT (DARK) will be colder (warmer) than the control, leading to changes in air-sea heat fluxes and potentially altering ventilation and deepwater formation rates. We deliberately use relatively large perturbations to the reflectance parameter in order to create changes in surface air temperature over Antarctica that are clearly distinguishable from noise.
The sensitivity experiments are started from the same initial conditions as the original 1850 control simulation, and the simulations extend 500 years using 1850 atmospheric composition.
We show in Figure  with that cooling extending at depth well into the Northern Hemisphere. This is consistent with both an enhanced rate of AABW formation, and AABW that is colder than in the control simulation. The reduced albedo case (Figure 24d) shows changes of opposite sign, but smaller in both amplitude and spatial extent. The reduced ventilation and AABW formation indicate that less heat is released from deeper layers of the SO, leading to subsurface warming. There is a reduction of the (cold) AABW penetration into the rest of the world ocean.
The middle row of Figure 24 shows changes in ocean salinity. The increased albedo case (Figure 24b) shows increased salinity in the near-surface layers around the coast of Antarctica. This is associated with an increased rate of sea ice formation and brine rejection associated with the cooling and is also consistent with increased vertical mixing. There is not a strong global scale signal of salinity change. For the reduced albedo case (Fig. 24e) there is a stronger signal of reduced salinity in the near-surface and increased salinity at depth. The reduced ventilation and AABW formation indicate that the warmer, saltier water at depth is not effectively mixed with the colder, fresher waters in the near-surface layers.
The bottom row of Figure 24 shows changes in zonal mean density. For the increased albedo case (Figure 24c) there is a very clear signal of increase in density for the entire SO, and the deepest layers of the rest of the world ocean. This is consistent with enhanced AABW formation, and AABW that is colder than in the control simulation. Comparing the temperature changes to the salinity changes, it is clear that most of the increase in density is a result of colder water rather than more saline water, with the exception of the near-surface layers off the coast of Antarctica. There is also an increase in density in the thermocline extending from approximately 50 o S to 50 o N. This may be related to colder Antarctic Intermediate Water (AIW) than in the control simulation, associated with the enhanced airsea heat fluxes in the higher latitudes of the SO.
For the reduced albedo case (Figure 24f), there is a reduction of density in the upper ocean at high latitudes of the SO, associated largely with reduced salinity. The reduced rate of AABW formation reduces the mixing of more saline deeper ocean water with the fresher upper ocean, thereby effectively reducing salinity and density. There is also a reduction of salinity in the thermocline, extending from approximately 50 o S to 50 o N. This appears to be linked to fresher AIW in the reduced albedo simulation than in the control. Since the nearsurface waters in the SO are the source waters for AIW, a freshening in the SO, associated with reduced albedo over Antarctica, is then communicated to the global thermocline via altered characteristics of AIW.  (Böhm et al., 2015).
An important feature of the global ocean circulation is the Antarctic Circumpolar Current (ACC). This zonal current is related to the meridional gradient of density between the coast of Antarctica and the coast of South America. A stronger (weaker) gradient, with denser water along the coast of Antarctica, is associated with a stronger (weaker) zonal current.
Since the albedo change over Antarctica alters this meridional density gradient (see Figure   24, bottom panels), we would expect an impact on the ACC. Figure 27a shows the time series of the ACC, calculated as the annual-mean volume transport through the Drake Passage. The three curves differ sharply with the strongest (weakest) ACC associated with the BRIGHT (DARK) case. The variation is substantial, from approximately 110 Sv (1 Sv = 10 6 m 3 s -1 ) in the reduced albedo case to almost 180 Sv in the increased albedo case. This substantial sensitivity points to the important influence of the energy balance over Antarctica on the ACC. An intriguing feature is the marked near-centennial scale variability in the control simulation that appears weaker in the perturbed simulations. This variability has resemblance to variability seen in other versions of GFDL models, and that has been hypothesized to play a role in the interpretation of recent decadal trends in sea ice and SST in the SO  .
As we have seen from zonal mean temperature and age, the responses to the albedo changes over Antarctica propagate through the world ocean. We show in Figure 27b the time series of the global mean net radiative balance at the top of the atmosphere. There is a clear dependence on Antarctic albedo, with the reduced (increased) albedo case having the most (least) positive imbalance at the top of the atmosphere. These differences persist throughout the 500 years of the simulations. Since the world ocean is the largest reservoir of heat in the climate system, the net radiative imbalance at the top of the atmosphere is usually linked with changes in the global mean, volume mean ocean temperature. We show time series of global mean, volume mean ocean temperature in Figure 27c. There is a wide range in long-term ocean heating or cooling, dependent on Antarctic albedo. From one perspective this is not surprising since Antarctic albedo contributes to global albedo and hence to global net radiative balance. However, it is not simply the direct impact of local albedo changes over Antarctica that is crucial; the related changes in ocean heat uptake, associated with fluxdriven changes in SO deepwater formation, play a fundamental role in the long-term changes in global mean ocean temperature.

Summary and Discussion
We have documented the development and simulation characteristics of a newly developed coupled ocean-atmosphere-land-sea ice model (Seamless system for Prediction and EArth system Research, "SPEAR"). This model is intended for studies of seasonal to multi-decadal climate variability, predictability and change, as well as experimental real-time seasonal to decadal prediction. This new modeling system incorporates many recent developments at GFDL, including improvements in the atmosphere, ocean, land and sea ice components, and the capability for variable resolution grids. The biases and simulation characteristics of the SPEAR model are significantly improved relative to corresponding characteristics in the previous generation GFDL seasonal to decadal prediction models: CM2.1  and FLOR (Vecchi et al., 2014). We also present the design characteristics and simulation characteristics of new 30-member ensembles of climate change simulations using SPEAR_LO and SPEAR_MED, extending from 1851 to 2100 using observed and projected changes in atmospheric composition and radiative forcing.
We also examined the sensitivity of aspects of the Southern Ocean (SO) and global climate system to the surface energy balance over Antarctica. Part of the motivation for these sensitivity experiments came from the early stages of the SPEAR development process during which the model produced little or no deepwater in the Southern Ocean (SO). This appeared to be related to relative warmth and freshness in the upper ocean around Antarctica and led to relatively strong stratification in the SO. This deficiency contributed to substantial long-term drifts in the world ocean; the reduced ventilation led to a buildup of heat in the interior of the SO from the transport of warmer waters from the Atlantic into the SO.
Multiple sensitivity tests were conducted to explore this bias involving ocean mixing and surface energy balance terms. We found that changing the heat balance over Antarctica by altering aspects of the snow albedo formulation had a strong impact on SO deepwater formation. We explored this in more detail by conducting additional sensitivity tests in which we altered terms in the snow albedo formulation. When albedo is increased surface air temperature is reduced over Antarctica relative to a control simulation. The cooled air flows out over the Southern Ocean (SO) and, due to an enhanced ocean-atmosphere temperature gradient, extracts more heat from the ocean than in a comparable control simulation. This cools the near-surface ocean, increasing upper ocean density through both cooling and enhanced sea ice formation and brine rejection, increases ventilation and deepwater formation, and leads to the formation of colder, denser Antarctic Bottom Water (AABW) that spreads at depth through the world ocean. More locally in the SO, the increased ocean density near the coast of Antarctica increases the meridional density gradient between Antarctica and South America, increasing the strength of the Antarctic Circumpolar Current (ACC) as measured by volume flow through the Drake Passage. The deep cooling of the world ocean leads to a substantial reduction in global ocean heat uptake in long control simulations, as well as reduced overall model drift. Generally opposite changes occur for the case of reduced albedo ("DARK"), although the amplitude of the impacts is somewhat smaller.
One of the additional intriguing findings in this study is that this model is able to simulate dense water formation through shelf processes that lead to dense water flowing down along the continental slope into the deepest layers of the SO (see Supporting Information Figure S10). This process occurs in addition to open ocean convection in the SO.
Formation of AABW from shelf waters has been a challenge for coarse-resolution models, and we find this development to be encouraging. This may be the result of improved numerical algorithms and/or the use of hybrid coordinates in MOM6 and remains the subject of ongoing exploration. The SO is difficult to model due to the importance of small-scale processes, such as mesoscale eddies, small-scale mixing processes, topographic variations, and the difficulty in simulating clouds and the surface radiation balance over the SO and Antarctica.
This study does not include the effects of freshwater from melting of Antarctic land ice (e.g., ice shelf, ice sheet) that could contribute to surface freshening and increased stratification. Inclusion of that process (Bronselaer et al., 2018) may well lead to reduced AABW formation but is beyond the scope of the present study.
The SPEAR model will be the primary tool for GFDL participation in a number of national and international prediction and attribution activities over the next several years.
These include the North American MultiModel Ensemble (NMME) for seasonal prediction (Kirtman et al., 2014) which uses models from five modeling centers across North America ©2020 American Geophysical Union. All rights reserved.
(https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/north-american-multimodel-ensemble). Real-time seasonal predictions are conducted once per month. SPEAR will also be used for decadal prediction activities contributing to CMIP6 and with ongoing predictions being made available through the UK Met Office as a WMO Lead Centre for Interannual to Decadal Prediction (https://www.metoffice.gov.uk/research/climate/seasonalto-decadal/long-range/wmolc-adcp). In addition, ensembles of SPEAR at high resolution are being used for the prediction and projection of changes in climate extremes.
The current version of the atmospheric component of SPEAR has 33 vertical levels.
Work is underway exploring the simulation characteristics of a version of this model with 65 levels. This version has a substantially improved simulation of the climatological characteristics of sudden stratospheric warmings and their connection to near-surface climate, but the computational cost of the model is also doubled. This version of SPEAR will be used for exploration of the impact of the stratosphere on seasonal to decadal climate variability and predictability.        (Cavalieri et al., 1996(Cavalieri et al., , climatology for 1990(Cavalieri et al., -2017, red is SPEAR_LO and blue is SPEAR_MED.      Table 1 for FLOR reference), (c) calculated from years 11-90 of SPEAR_LO, and (d) calculated from years 11-90 of SPEAR_MED. The green ellipses denote the centers of the teleconnection patterns from the reanalysis in (a).                 Vecchi et al. (2014)