Projections of physical conditions in the Gulf of Maine in 2050

, The Gulf of Maine (GoM) is currently experiencing its warmest period in the instrumental record. Two high-resolution numerical ocean models were used to downscale global climate projections to produce four estimates of ocean physical properties in the GoM in 2 0 5 0 for the “business as usual” carbon emission scenario. All simulations project increases in the GoM mean sea surface temperature (of 1.1 (cid:2) C–2.4 (cid:2) C) and bottom temperature (of 1.5 (cid:2) C–2.1 (cid:2) C). In terms of mean vertical structure, all simulations project temperature increases throughout the water column (surface-to-bottom changes of 0 .2 (cid:2) C– 0 .5 (cid:2) C). The GoM volume-averaged changes in temperature range from 1.5 (cid:2) C to 2.3 (cid:2) C. Translated to rates, the sea surface temperature projections are all greater than the observed 1 00 -year rate, with two projections below and two above the observed 1982–2 0 13 rate. Sea surface salinity changes are more variable, with three of four simulations projecting decreases. Bottom salinity changes vary spatially and between projections, with three simulations projecting varying increases in deeper waters but decreases in shallower zones and one simulation projecting a salinity increase in all bottom waters. In terms of mean vertical structure, salinity structure varies, with two simulations projecting surface decreases that switch sign with depth and two projecting increases throughout the (subsurface) water column. Three simulations show a difference between coastal and deeper waters whereby the coastal zone is projected to be systematically fresher than deeper waters, by as much as 0 .2 g kg –1 . Stratification, 5 0 m to surface, is projected to increase in all simulations, with rates ranging from 0 . 00 3 to 0 . 00 6 kg m –4 century –1 which are lower than the observed change on the Scotian Shelf.The results from these simulations can be used to assess potential acidification and ecosystem changes in the GoM.


Introduction
The Gulf of Maine (GoM) is currently experiencing one of the fastest rates of warming of any ocean ecosystem . The GoM 2050 International Symposium held November 4-8, 2019, brought together members from the scientific, fishing, environmental, business, and other communities to discuss climate change and its potential impacts on the GoM. Presentations at the Symposium focused on projected changes to the physical environment and the related effects on ocean ecology and acidification, sea level, coastal communities, and commerce. High-resolution numerical ocean circulation models provided the climate projections for the physical environment that were used in the analyses of possible changes to the ocean's ecology and acidification level. This article describes the details of how these model simulations were performed and their key projections of changes to the physical environment of the GoM.
To understand possible future climate changes in the GoM, placing the region in the context of the large-scale ocean circulation is important. The Newfoundland/Labrador Shelf, Gulf of St. Lawrence, Scotian Shelf, and GoM form an interconnected system along the eastern seaboard of Canada and the United States. The circulation in the region is characterized by a general northeast-southwest flow of water from the Labrador and Newfoundland Shelf areas through the Gulf of St. Lawrence, Scotian Shelf, and GoM to the Mid-Atlantic Bight (Figure 1; see, e.g., Xue et al., 2000;Hannah et al., 2001;Han et al., 2008;Wu et al., 2012;Brickman et al., 2018). The region off the shelf is the confluence zone between the warm northeastwardflowing Gulf Stream and the cold southwestwardflowing Labrador Current (Loder et al., 1998). These two currents interact at the tail of the Grand Banks (south of Newfoundland) resulting in subsurface east-to-west flows along the shelfbreak that affect ocean variability downstream on the Scotian Shelf and GoM. Ocean properties in the GoM are affected by the Nova Scotia (NS) Current, which brings colder, fresher waters from the Gulf of St. Lawrence, and which merges into the GoM Coastal Current, creating the general counterclockwise circulation found in the region. Another important inflow is the warm/salty subsurface inflow through the eastern Northeast Channel (NEC;Smith, 1983). These inflows are balanced by outflows at the Great South Channel and the western part of the NEC (along the eastern flank of Georges Bank). Other factors influencing ocean properties in the GoM are warm core rings (from the Gulf Stream) and local effects like river inputs and interaction with the atmosphere.
The characteristics of the circulation and surface inputs result in characteristic spatial patterns and variability in the GoM and regions upstream. For example, the bottom conditions change from cold/fresh on the eastern Scotian Shelf and points east to warm/salty from the central Scotian Shelf and west into the GoM (Figure 1). This change is due to the influence of the off-shelf warm/salty recirculation waters that penetrate via deep channels onto the shelf regions, with the NEC being the primary conduit into the GoM. Brickman et al. (2018) studied the variability of these deep waters and found that this was controlled by the interaction of the Gulf Stream and Labrador Current at the Tail of the Grand Banks that results in east-to-west propagating eddies. During the last 10-15 years, this process has resulted in a dominance of warm/salty Gulf Stream eddies that are the major contributor to the recent warming trends in the region's deep waters Hebert et al., 2018;Galbraith et al., 2019;Friedland et al., 2020).
While the temperature (T), circulation, and other physical aspects of the GoM are changing on decadal and longer time scales (Shearman and Lentz, 2010;Pershing et al., 2015;Kavanaugh et al., 2017;Brickman et al., 2018;Chen et al., 2020), the GoM also experiences shorter term fluctuations including marine heat waves (Mills et al., 2013;Scannell et al., 2016;Pershing et al., 2018). Marine heat waves can be defined as periods when daily sea surface temperatures (SSTs) are above the seasonally varying 90th percentile (threshold) for at least 5 consecutive days (Hobday et al., 2016). The GoM has experienced record SSTs during the last decade. For example, during the 2012 heat wave, the SST was 1.3 C warmer than the long-term average and on par with some end-of-century climate projections (Mills et al., 2013). Observational analyses and model simulations indicate that the 2012 SST warming event was driven primarily by the atmosphere through anomalous air-sea heat fluxes (Chen et al., 2014. Relative to today's climate, marine heat waves are projected to increase in the future primarily through an overall warming of the ocean, that is, through a mean shift  in the climate, rather than due to changes in the variability of temperature anomalies Frölicher et al., 2018;Oliver, 2019). The Intergovernmental Panel on Climate Change (IPCC) reports on the future climate projections from an ensemble of coupled atmosphere-ocean earth system models (ESMs). The computational demands of these simulations require that the ocean model component of these ESMs have horizontal resolutions O (100 km), which are considered to be low-resolution ocean models. Indeed, the processes affecting the GoM (described above) could not be simulated accurately by these models. To properly resolve the scales of interest in our region requires a technique called dynamic downscaling. Briefly, dynamic downscaling uses future climate output from the IPCC coarser resolution global ESMs to drive higher resolution regional models of the ocean, thus producing high-resolution future climate projections. This article focuses on the outputs from two techniques for downscaling future climate simulations to regions of the North Atlantic Ocean-outputs which provided the basis for other analyses associated with the GoM 2050 Symposium and reported in this special feature.

Methods
The main tools for projecting future climate conditions are coupled global atmosphere-ocean-sea ice-land models and ESMs that also include representations of ecosystems and chemistry, which are run by numerous international institutes. The IPCC coordinates the simulation, analyses, and reporting of these future climate simulations in a series of "Coupled Model Intercomparison Projects" (CMIP) of which CMIP5 (the fifth CMIP) forms the basis of the most recent IPCC report (Taylor et al., 2012). (The sixth CMIP has been completed with report due in 2021.) Future projections for the CMIP5 simulations used Representative Concentration Pathway (RCP) emission scenarios, which specify different time trajectories for concentrations of greenhouse gases based on socioeconomic factors (IPCC, 2014). Four RCP scenarios were developed for CMIP5 (2.6, 4.5, 6.0, and 8.5), where the number refers to the radiative forcing in Wm -2 in 2100. Two RCPs (4.5 and 8.5) were used by the models reported herein, although we focus primarily on simulations that used the RCP8.5 scenario. RCP4.5 is an intermediate scenario in which emissions peak around 2040, then decline until 2100. RCP8.5, often referred to as "the business as usual" scenario, is the most extreme case in which the radiative forcing increases through the 21st century.
The horizontal resolution of the CMIP5 ocean models is considered to be coarse, ranging from about 0.5 -2 (in longitude) which translates into a resolution of 50 km or greater in the northwest Atlantic Ocean. This coarseness presents a challenge in this region as the main current systems (e.g., the Gulf Stream and Labrador Currents) are not properly resolved and thus not accurately simulated in these models. In addition, CMIP5 models do not resolve narrow coastal currents, ocean eddies, and interactions with the complex bathymetry along the northeast U.S. shelf. The result is large biases (errors) in these model simulations of the present-day climate, which reduces confidence in the future climate simulations. More precisely, the position of the Gulf Stream is typically too far north, which results in a warm bias in the GoM and off-shelf region.
One of the implicit assumptions of future climate modeling is that while models may have weaknesses in simulating aspects of the present climate, the differences between their future climate projections and present climate simulations contain useful information. This assumption is the basis of the "delta method." The first step in the delta method is to select a common period from the simulations and in the real world (e.g., 1976-2005). Then the difference, or delta, between this period and the target period (in this case, 2050) is computed for variables of interest in the simulations. The delta values can be used directly, to portray projected changes and their significance, or can be added to the observed values for the reference period from the real world, if absolute values of variables are required. This approach has been used to support several ecosystem projections for the GoM (e.g., Hare et al., 2010;Kleisner et al., 2017;Le Bris et al., 2018;Greenan et al., 2019).
A finer scale representation of the present-day climate and how it changes in the future can be obtained using the dynamic downscaling method, in which output from the low-resolution CMIP models is used to force a highresolution model for a portion of the globe. Two modeling groups have applied this dynamical downscaling approach over domains that include the GoM. The National Oceanic and Atmospheric Administration's (NOAA) Physical Sciences Laboratory used output from three different CMIP5 models-Geophysical Fluid Dynamics Laboratory (GFDL), Institute Pierre Simon Laplace (IPSL), and Hadley Global Environment Model (HadGEM)-to drive a highresolution ocean model that extended from the Gulf of Mexico to Newfoundland. Using the three different global models provides a way of capturing some of the range of possible future conditions. We will refer to the output from this model as the Regional Ocean Modeling System (ROMS) simulations. Canada's Department of Fisheries and Oceans used a similar procedure for a highresolution ocean model of the North Atlantic Ocean. This regional model was forced with the average output from six IPCC models run under both the RCP4.5 and 8.5 scenarios. We will refer to these as the Bedford Institute of Oceanography North Atlantic Model (BNAM) simulations. Technical details on both the BNAM and ROMS approaches are described below.

BNAM methods
The Fisheries and Oceans Canada future climate modeling was done using BNAM (Brickman et al., 2016Wang et al., 2019). BNAM is a high-resolution (1/12-degree) model of the North Atlantic Ocean, based on the NEMO-OPA code (Madec, 2008), on a domain of 7-75 N latitude and 100 W-25 E longitude ( Figure S1). The z-level model has 50 vertical levels, partial cells for the bottom layer, and a horizontal resolution in the GoM region of about 6 km. BNAM was designed specifically to resolve the interaction of the Gulf Stream and Labrador Current at the tail of the Grand Banks, which is required to simulate changes in Maritime Canadian waters, including the GoM. BNAM is run in present and future climate modes. Forcings for the BNAM future climate simulations were derived as anomalies from the mean of an ensemble of six IPCC coupled atmosphere-ocean future climate runs for two future periods, 2055 (2046-2065) and 2075 (2066-2085), and two RCPs, 8.5 and 4.5. The four future climate scenarios also included projections of future river runoff and a representation of the expected increase in melting of the Greenland glacier (the latter has a noticeable effect).
The present ocean climate was simulated using the (repeat cycle) Co-ordinated Ocean-Ice Reference Experiments (CORE) Normal Year atmospheric forcing (Large and Yeager, 2004). Future climate anomalies were added to the present climate to create the four future climate forcings. The resulting ocean model simulations produced climatologies for the future periods 2055 and 2075 for RCPs 8.5 and 4.5. Results are typically presented as spatial fields of future climate changes, that is, as differences between the model future and present climates. This approach allows the model output to be added to present climate fields derived from a variety of sources (Delta method).

ROMS methods
The effects of climate change on the northwest Atlantic were investigated using ROMS (Shchepetkin andMcWilliams, 2003, 2005). The version used here, configured by Kang and Curchitser (2013), has a horizontal grid spacing of 7 km and 40 vertical levels. The model domain covers the Gulf of Mexico and extends along the entire U.S. east coast from Florida to Newfoundland, including the entire GoM ( Figure S1).
A control (CTRL) ROMS simulation was performed using observationally based fields at the surface and along the side boundaries in the ocean, with radiation conditions for flow out of the domain and nudging of temperature, salinity (S), and in-flowing currents as a function of depth at the boundary (see Alexander et al., 2020, for more details). The initial conditions and oceanic boundary forcing were derived from 5-day averages from the Simple Ocean Data Assimilation (SODA v2.1.6; Carton and Giese, 2008), 6-h surface forcing from CORE version 2 (Large and Yeager, 2009), and daily freshwater flux from rivers from the continental discharge data set (Dai et al., 2009). The CTRL simulation represents the observed climate over the period 1976-2005 and well simulates the mean path of the Gulf Stream Curchitser, 2013, 2015), the circulation in the GoM , and temperatures on the Northeast U.S. Shelf .
The ROMS climate change simulations were also performed using the delta method. The delta values used to derive the initial and boundary conditions in the future ROMS simulations were obtained from three CMIP5 models: GFDL ESM, GFLD ESM2M; IPSL climate model, CM5A-MR, and HadGEM climate configuration. We will refer to these as the GFDL, IPSL, and HadGEM simulations (or simply GFDL, IPSL, and HadGEM). The three ESMs differ in terms of how they simulated the present-day climate (e.g., strength of the Gulf Stream) and how strongly they responded to greenhouse gases. The deltas were computed in each of the three ESMs by subtracting the long-term monthly mean values during the period 1976-2005 from those during 2070-2099 and then adding them to the CTRL in order to perform three independent future simulations. The ROMS simulations do include the estimates of river runoff and how they change in the future but not glacial melt. Additional information about the ROMS CTRL and future simulations, including their representation of climate change in the northwest Atlantic Ocean, can be found in Alexander et al. (2020) and Shin and Alexander (2020).
Because we used the 2070-2099 period for the future climate forcing, the ROMS experiments needed to be adjusted to estimate the changes in 2050. We explored several ways to accomplish this adjustment and tested them using the CMIP5 archive, which provided continuous monthly time series over the 21st century, including midcentury values. One method that worked well for all variables was to use the global average change in radiative forcing. For the RCP8.5 scenario, the changes in forcing from the present climate  to 2050 and 2085 are 3.1 and 5.68 W m -2 , respectively. So, the differences between the CTRL and climate change simulations presented here were scaled by 0.546 (3.1/5.68) to estimate the changes that occurred by 2050.

Results
Together, the BNAM and ROMS simulations provide five different views of the future state of the GoM in 2050 (RCPs 8.5 and 4.5 from BNAM, plus three for RCP8.5 from ROMS). These are expected to differ due to the sensitivity of their respective global models to carbon dioxide levels and, in the BNAM simulations, to differences in carbon emission pathways. Some of the differences may also be due to the two distinct modeling approaches. For example, the BNAM simulations include a detailed treatment of the projected Greenland glacier melt which does show a downstream effect along the eastern seaboard of Canada and the United States.
In this section, we present outputs from the BNAM and ROMS simulations for key ocean variables. We focus on providing results that summarize the horizontal spatial patterns (represented by annual average changes in stratification and surface and bottom T and S) and their variability in space (across the GoM) and time (month). The vertical dimension is represented by vertical profiles of T and S averaged over the GoM and their projected changes in 2050. We provide a brief summary of the models' projected changes in circulation in the region. Below, we present only the four simulations that used the RCP8.5 "business as usual" scenario. The patterns in the BNAM RCP4.5 and RCP8.5 simulations are similar, differing mainly in the intensity of the response. Thus, presenting the details of the RCP4.5 simulation is not highly informative. Instead, we utilize the BNAM RCP4.5 results in the last section to provide comparison from a more optimistic emissions scenario.

SST
The four model simulations that run under the RCP8.5 scenario all show increases in annual averaged SST throughout the GoM with spatial variability of about 0.5 C ( Figure 2). The BNAM simulation suggests warming of 1 C-1.5 C above its present climate baseline. The three ROMS simulations suggest warming of the annual mean SST over the entire GoM, ranging from 1.5 C to 2.0 C, 1.75 C to 2.25 C, and 2.25 C to 2.75 C in the simulations driven by the GFDL, IPSL, and HadGEM climate models, respectively ( Figure 2). This ordering is consistent with the overall sensitivity of the three climate models: The Hadley model has the strongest global warming and warmest temperatures over North America and the Arctic, followed by IPSL and then GFDL. The projections for SST change can be averaged over the GoM, converted to a rate of change, and compared to estimates based on data. Fernandez et al. (2015) and Pershing et al. (2015) report changes for two time periods: long-term, 1900-2013 ¼ 0. 6 C century -1 ; and recent, 1982-2013 ¼ 3.0 C century -1 . The GoM average change for BNAM is 1.1 C, which translates to a rate of 1.8 C century -1 (using the central value of the model's present climate period of 1995). For the ROMS projections (present climate central value of 1990), we get 2.7 C century -1 , 3.2 C century -1 , and 3.8 C century -1 for GFDL, IPSL, and HadGEM, respectively. Thus, the BNAM and GFDL projections lie between the long-term and recent climate trends, while both IPSL and HadGEM exceed the recent trend.

Sea surface salinity (SSS)
Changes in SSS (Figure 2) provide an example of how differences in ocean model future climate forcing and processes can result in variability in intermodel projections. The ROMS simulations exhibit substantial differences among the three simulations, with a decrease in the GFDL simulation and HadGEM but an increase in IPSL. The freshening is especially strong in GFDL, resulting from a strong decrease in salinity over the whole North Atlantic and stronger southeastward coastal currents on the Scotian Shelf, both of which lead to fresher water being transported to the GoM Shin and Alexander, 2020). The BNAM simulation indicates the strongest salinity changes with a relatively uniform spatial decrease in the range of -0.2 to -0.3. (NB: Following convention, we report model salinity, itself dimensionless, without units.) The stronger BNAM SSS changes are likely due to the increase in river input to the GoM in 2050 (5.6% increase in freshwater transport from rivers), precipitation, as well as a long-range effect due to the southward advection of freshwater from the increased Greenland glacier melt included in that model.

SST and SSS spatial variability as a function of season
To portray the variability of SST and SSS as a function of space and time, box-and-whisker plots were created from the monthly individual grid points in the GoM region ( Figure  3), where the ("entire") GoM was defined  as the [longitude, latitude] box [70.875-65.375 W, 40.375-45.125 N] (we excluded depths > 350 m). For SST (Figure 3a), the ROMS model warming varies with the seasons, with the strongest warming occurring in late summer- early fall for the IPSL and GFDL simulations, enhancing the seasonal cycle of SST, in agreement with the analysis of CMIP5 models by Alexander et al. (2018). The enhanced warming occurs when the mixed layer is shallow in summer, and thus, the heating from the atmosphere resides closer to the surface. Surface temperature changes in BNAM exhibit a minimum during the summer season, with greatest spatial variability from late summer to early winter. For SSS (Figure3b), the ROMS annual cycle of surface salinity response is opposite to that of its temperature, with the strongest decreases occurring in summer in all three simulations, when even the IPSL values are negative. BNAM surface salinity changes (the freshest) exhibit the same seasonal variability, with greatest spatial variability from late summer to early winter. Analysis of the spatial patterns showed that this variability is due to a shift in the NS current (NSC) as it circulates around the western tip of NS (see Figure 1). This current transports the spring freshwater pulse from the Gulf of St. Lawrence into the GoM arriving late in the year. In the future climate projection, the position of the freshwater input shifts away from the coastline toward the west, which manifests itself as large spatial variability when differenced from the present climate pattern.

Bottom temperature and salinity
All of the simulations show increased temperature and salinity in the deep basins of the GoM (Figure 4). BNAM-projected changes in bottom temperatures range from 1.5 C to 3 C in the deeper waters of the GoM and from 0.5 C to 1.5 C in the shallower coastal waters and on Georges Bank (Figure 4). Changes in bottom salinity show an interesting pattern with freshening in the shallower coastal waters and on Georges Bank by about -0.25, while the deeper waters show an increase in salinity of about 0.25. This pattern indicates that the deeper waters are influenced by intrusions of warm salty off-shelf waters (through the NEC), while the shallower waters are affected by BNAM's projection of increased river inputs and precipitation in the GoM region and advection of freshwater from the Arctic.
ROMS-simulated bottom temperature changes indicate warming on the order of 1.0 C-2.0 C for GFDL and IPSL and 2.0 C-2.75 C for HadGEM by 2050 (Figure 4). The warming varies with depth among the three models. The warming is greatest at depth in GFDL and enters the GoM via the deep NEC. The maximum warming occurs higher in the water column and intersects the bottom on the shallower parts of the shelf around the Gulf in HadGEM and to a lesser degree in IPSL.
The bottom salinity in both the present day and the future is strongly influenced by the complex bathymetry in the GoM. In the present day, relatively salty water enters the Gulf via the NEC and is ringed by less saline water in shallower portions of the Gulf ( Figure S14). A similar process occurs to varying degrees in the climate change simulations. The increase in salinity is confined to the vicinity of the NEC and the adjacent Gordons and Crowell basins in GFDL but is much more expansive in the other two simulations. The shallower parts of the GoM exhibit a decrease in salinity in GFDL (similar to the BNAM result) and to a much lesser degree in HadGEM, in concordance with the freshening at the surface in these two models.

Vertical profiles of temperature and salinity
Spatial mean, annually averaged, vertical profiles of temperature and salinity were computed for the top 200 m of To see whether the coastal zone behaved differently, the spatial mean, annually averaged, vertical profiles of temperature and salinity were computed for the waters inside the 100-m isobath, exclusive of Georges Bank (Figure 5c and f). For BNAM, the coastal zone profiles are essentially the same as those for the top 100 m of the basin (Figure 5a and d), indicating the degree to which the off-shelf intrusions which enter through the NEC flood the entire deeper layers of the GoM in that model. For ROMS, the profiles show a clear freshening in the coastal zone (below 10 m) relative to the regional mean, with a slight tendency toward cooler waters. The freshening is strongest in the GFDL simulation (*0.1) with the IPSL and HadGEM simulations showing similar but weaker (*0.05) tendencies. For the latter two simulations, this change toward relatively fresher waters in the coastal zone is not strong enough to erase the overall increase in salinity projected to occur throughout the water column for those simulations.

Changes in stratification
There are numerous definitions of vertical stratification. Here, we use the commonly used shelf definition (Li et al., 2015;Hebert et al., 2018) as the difference in density between the surface and 50 m, in kg m -3 , expressed as a positive number (i.e., r(50)-r(0), where r is the density). Figure 6 shows the spatial plots of the annually aver-  With respect to the seasonality of the stratification changes (Figure 7), all models show a seasonal cycle with a peak around months 7-9. This seasonal cycle is consistent with the models' present climate seasonal cycle of stratification (not shown) and other sources (see, e.g., Li et al., 2015, their Figure 5;Johnson et al., 2018, their Figure 6). With respect to the amplitude of the changes, the rankings reflect the mean spatial distributions: GFDL, BNAM, HadGEM, and IPSL, from lowest to highest.
To put these changes into perspective relative to trends from regional data, the projected changes in 2050 can be converted to a rate of change based on the central value of the models' present climate period (as done above). Using the spatial means for stratification (and dividing by the 50-m depth range to convert to vertical gradient), we project rates of change of approximately 0.003, 0.004, 0.005, and 0.006 kg m -4 century -1 for GFDL, BNAM, HadGEM, and IPSL, respectively. The 70-year trend from Hebert et al. (2018) for the Scotian Shelf region is 0.014 kg m -4 century -1 . Thus, the models are predicting a continuing increase in stratification in the region, with rates similar to but less than the present rate reported for the Scotian Shelf.

Changes in circulation
Quantifying and summarizing changes in velocity, and circulation, is less straightforward than for scalar variables like temperature and salinity. All three ROMS simulations indicate a weakening of the Atlantic meridional overturning circulation, with related reduction (*25%) in the strength of the Gulf Stream . Although the simulations project differences in the Gulf Stream position, these differences did not directly affect the circulation in the GoM. The ROMS projections do indicate an increase in the counterclockwise circulation in the region, but this is due to increases in flow entering the GoM from the east via the NEC and NS current (NSC). By contrast, the BNAM simulations project reduction in this circulation due to decreases in transport into the GoM through the NEC and NSC, with the fraction of transport into the GoM from the NSC (NSC/[NEC þ NSC]) decreasing from 0.55 to 0.46 by 2050 (see Brickman et al., 2016).

Discussion and summary
This article has reported on four simulations of the physical environment of the GoM in 2050, based on the business-as-usual RCP8.5 scenario. These simulations were computed by two downscaling ocean modeling systems: ROMS and BNAM. The ROMS system performed three simulations based on output from the GFDL, IPSL, and HadGEM ESMs, which allowed for variability in its forcing. The BNAM system was forced by the average of six ESMs which serves to smooth out variability in the ESM forcing. The results (Figures 2-7) provide an estimate of the expected changes in the GoM plus its variability. We find that all four simulations project increases in surface and bottom temperatures (Figures 2 and 4) and, based on the spatial mean temperature profiles ( Figure  5), this statement can be extended throughout the water column. Calculations of the volume-averaged GoM temperature increases are 1.5 C for BNAM and GFDL, 1.7 C for IPSL, and 2.3 C for HadGEM. While these order in the same way as the mean SST changes, over the entire GoM, BNAM projects a volume average temperature 0.4 C higher than its SST value, while GFSL, IPSL, and HadGEM volume averages are lower by 0.1 C, 0.2 C, and 0.1 C, respectively, than their SST values. This finding is consistent with the mean T profiles shown in Figure 5.
How salinity may change is not so clear. For example, for SSS (Figure 2), the projections generally support a freshening, but one model (IPSL) projects an increase in salinity. Similarly for bottom salinity, the BNAM simulation projects increased salinity in the deeper waters but a freshening in shallower waters (<100 m), which is supported by GFDL and HadGEM to a lesser extent, but not IPSL. The mean vertical salinity profile from BNAM projects a freshening for the top 70 m and salinification below (for both the entire and coastal GoM), while the ROMS results are more variable: The GFDL simulation switches from negative to positive salinity changes at *120 m, while the HadGEM and IPSL simulations project increases in salinity throughout the water column. The ROMS results also show a difference between the coastal and deeper waters whereby the coastal zone is projected to be systematically fresher than deeper waters.
With respect to space-time variability (Figure 3), all model simulations project positive changes for SST throughout the year but with variability in the monthly cycle. For SSS, the model results suggest that same seasonality (with minima in late summer), but the sign of the change varies between simulations: BNAM and GFDL are uniformly fresher throughout the year, while HadGEM and IPSL change sign during the course of the year. Regarding circulation changes, the ROMS and BNAM simulations indicate opposite tendencies for the counterclockwise circulation in the Gulf, with BNAM projecting a weakening and ROMS projecting a strengthening of this circulation. (NB: Views of the ROMS velocity and scalar fields are available at https://psl.noaa.gov/ipcc/roms/.) Saba et al. (2016) also investigated the response of the GoM to climate change, using a global coupled modeling system with different horizontal resolutions. They found that the response varied with resolution, with the increase in temperature in the GoM strongest in their highest resolution 10-km ocean model. The strong warming in the GoM is consistent with the results from the ROMS and BNAM simulations shown here. However, Saba et al. (2016) attributed the enhanced warming in the GoM to a northward shift of the Gulf Stream, while in the ROMS simulations, this warming was attributed to circulation changes . The projections of climate change in the GoM differ. Why does this occur? Three main factors cause differences in climate projections: (1) the amount of greenhouse gases that will be released into the atmosphere, (2) differences in how models are formulated, and (3) variability that arises naturally in the climate system (Hawkins and Sutton, 2009). The uncertainty of the first has been tested using different climate change scenarios, and the second can be tested by using different climate models. However, natural variability can complicate isolating humaninduced climate change, especially for small regions of the globe such as the GoM. Natural variability can be seen in long-observed time series, with periods of both rapid increases but also decreases in temperature (Figure 8), with natural variability likely contributing to the very warm temperatures that occurred recently. In addition, changes over decades or longer due to natural variability were often assumed to be small. However, there are very long time scale fluctuations that occur due to slow changes in ocean currents, especially the meridional overturning circulation, which strongly influences the North Atlantic Ocean.
A method for examining natural variability in a changing climate has recently been developed in which a large set of climate model simulations are performed using the same scenario and model but with different initial conditions for each simulation. The differences at the start of the integrations cause changes in the temperature, winds, currents, and other factors to vary, resulting in each simulation having a different evolution of the atmosphere, ocean, and sea ice. The average of all the simulations gives the model's response to climate change, and the spread among the simulations indicates the influence of natural variability. Results from these large ensembles of simulations are eye opening: Even over 50-year periods, natural variability could strongly influence regional trends, especially for variables such as precipitation and sea-level pressure (Deser et al., 2012(Deser et al., , 2014(Deser et al., , 2020Kay et al., 2015). Results from a large ensemble for SSTs in the GoM suggest inherent variability of about +0.5 C.
To place our projections in context, we present the mean change in GoM SST relative to historical conditions ( Figure 8). Because each of the simulations represents the mean climate state (*30-year average) around 2050, comparing the projections with the 30year average of the observations (black line) is most appropriate. We see the most recent 30-year period (black cross centered on 2002 in Figure 8), which included the very warm conditions of the past decade (2010-2018), contrasted against the cooler conditions of the 1990s. Note as well that the current "climate" of the GoM is only now approaching the mean conditions around 1950. Although the GoM has already experienced brief conditions like those indicated by the climate models, each of the downscaled projections clearly suggests a climate that is significantly warmer than what the GoM on average has experienced.
From Figures 2 and 8  Years like 2012 and 2016 that now seem remarkable would feel like merely warm years. The very slight 0.3 C difference between BNAM RCP4.5 and RCP8.5 underscores the insidious delay between when carbon is released into the atmosphere and when its effects are fully felt in the climate system. 3. ROMS GFDL: This is the coolest of the ROMS simulations, but it is still much warmer than the BNAM simulations. The 1.55 C change in mean temperature is only slightly below the temperature anomaly in 2016 (1.82 C). In this climate, 2012 and 2016 would be ordinary years, and extremely warm years would have temperatures well beyond historical experiences. 4. ROMS IPSL: The 1.89 C change in mean temperature corresponds almost exactly to the anomaly during the record year of 2012 (1.88 C). We would expect that this climate could produce years with a mean temperature anomaly more than 2 C above the baseline. 5. ROMS HadGEM: This simulation indicates the most extreme warming. The 30-year average would be 2.38 C above the baseline, almost 0.5 C above the warmest year recorded. In this climate, a cold year would feel like 2012.
The projections presented in this study offer a view of the potential future for the GoM. We cannot determine whether one of these is more likely than the others. However, all of these model projections show a GoM that is warmer than the recent 30-year climate, so this projection is likely robust. Time series, such as that shown in Figure  8, indicates a long-term upward trend but with significant interdecadal and interannual variability. Recent studies indicate that positive SST extremes are increasing in frequency over much of the global ocean, primarily due to shifts in the mean climate (Diffenbaugh et al., 2017;Alexander et al., 2018;Oliver, 2019). Thus, variations in the climate system in conjunction with the projected continuing warming trend will result in increasing extreme warm events in the GoM.

Data accessibility statement
The observational sea surface temperature data are available at https://www.ncdc.noaa.gov/oisst. Many of the plots from the Regional Ocean Modeling System simulations can be recreated using the NOAA Climate change web portal: https://www.psl.noaa.gov/ipcc/ roms/.
Model outputs used in the analyses presented here are provided in main figures and Supplemental Materials; additional model outputs are available upon request (contacts: DB, MAA, JDS, ZW).

Supplemental files
The supplemental files for this article can be found as follows: Figures S1-31. PDF