Changes in extreme regional sea level under global warming

An important contribution to future changes in regional sea level extremes is due to the changes in intrinsic ocean variability, in particular ocean eddies. Here, we study a scenario of future dynamic sea level (DSL) extremes using a strongly eddying version of the Parallel Ocean Program. This model is forced with atmospheric fluxes from a coupled climate model which has been integrated under the IPCC-SRES-A1B scenario over the period 2000-2100. Changes in 10-year return 5 time DSL extremes are very inhomogeneous over the globe and are related to changes in ocean currents and corresponding regional shifts in ocean eddy pathways. In this scenario, several regions in the North Atlantic experience an increase in mean DSL of up to 0.4 m over the period 2000-2100. DSL extremes with a 10-year return time increase up to 0.2 m with largest values in the northern and eastern Atlantic. 10


Introduction
From satellite measurements, it has been well established that global mean sea level has increased by about 3 mm/yr over the period 1993-2010 (Rhein et al., 2013;Church and White, 2011;Church et al., 2013;Watson et al., 2015). However, regional sea level trends are very inhomogeneous over the oceans, and range from about a 10 mm/yr increase in the western tropical Pacific to about 5 1986-2005), with major contributions provided by thermal expansion of ocean water and the mass loss of the major ice sheets and glaciers (Church et al., 2013;Slangen et al., 2014). For the highest radiative forcing scenario (RCP8.5), projected global sea level rise is between 0.52 and 0.98 m in 25 2100 (Church et al., 2013). Regional sea level changes projected for the North Atlantic show complex patterns that are partly caused by a weakening of the Atlantic Meridional Overturning Circulation (AMOC), by a shift in the path of the North Atlantic Current, and by changes in surface buoyancy fluxes (Landerer et al., 2007;Yin et al., 2010Yin et al., , 2009Pardaens et al., 2011;Kienert and Rahmstorf, 2012;Bouttes et al., 2013). Thermosteric sea level evolves with a pattern that reflects the reduced 30 heat transport to the North Atlantic due to changes in ocean currents (Yin et al., 2010;Pardaens et al., 2011). For example, DSL is rising near the North American continent because of a reduction in the AMOC causing a redistribution of ocean mass (Landerer et al., 2007;Yin et al., 2009Yin et al., , 2010Bouttes et al., 2013).
The spread in the projections of regional sea level change is largely determined by internal ocean 35 variability and model uncertainty. In Slangen et al. (2012) and Bordbar et al. (2015), the spread due to decadal-to-centennial variability is considered by looking at ensemble simulations using CMIP3 and CMIP5 climate models, respectively. It was shown that the (CMIP5) ensemble spread of the projected DSL is of the same order of magnitude as the globally averaged sea level rise (Bordbar et al., 2015). Several regions were identified where the forced sea level change signal is relatively 40 strong with respect to the internal variability, e.g. the Indo-Pacific part of the Southern Ocean and the eastern equatorial Pacific and hence may be detected earlier (Bordbar et al., 2015).
However, in all these model studies the strongest component of oceanic internal variability, i.e. that due to ocean meso-scale eddies, was not represented. Rectification processes due to eddies can lead to strong changes in mean ocean surface flows and their response to atmospheric forcing, in particular 45 in the Southern Ocean (Böning et al., 2008). In strongly eddying ocean models even new modes of low-frequency variability may appear, such as the multidecadal Southern Ocean Mode (Bars et al., 2016). Using the eddy-permitting (about 1/4 • horizontal resolution) version of the MIROC3.2 model, Suzuki et al. (2005) showed that representing ocean eddies provides a more detailed projection of regional sea level changes under the IPCC SRES-A1B scenario and that eddies are strongly involved 50 in regional sea level extremes. In addition, as demonstrated by Firing and Merrifield (2004) from observational data, a high background sea level superposed on the sea level change due to an arriving ocean eddy can lead to extreme local sea levels.
Eddies can also have a strong effect on the deep ocean circulation, in particular the AMOC. The study of Weijer et al. (2012) indicates that the AMOC in the strongly eddying (about 0.1 • horizontal 55 resolution) version of the POP model version is more sensitive to freshwater perturbations than the non-eddying version of the same model. Climate model studies on the projections of the AMOC with non-eddying ocean components show only an AMOC decline of 22% to 40% over the period 2000-2100, depending on the IPCC scenario (Weaver et al., 2013). Only two (out of 30) of these 2 Ocean Sci. Discuss., doi:10.5194/os-2016Discuss., doi:10.5194/os- -57, 2016 Manuscript under review for journal Ocean Sci. Published: 2 August 2016 c Author(s) 2016. CC-BY 3.0 License. models project a substantial decrease of the AMOC under the RCP8.5 scenario. However, in high-60 resolution ocean models strong variations in the AMOC strength lead to changes in ocean currents and eddy pathways, which induce an additional contribution to the variability in DSL and hence affect extreme DSL values (Brunnabend et al., 2014).
It is important to assess the role of eddies in projections of future regional sea level changes, in particular on the DSL extremes. In this paper, we study a scenario of future DSL change using the 65 strongly eddying version of POP as in Brunnabend et al. (2014), but now forced with atmospheric fields from a coupled climate model that evolved under the SRES-A1B scenario. We focus on the changes in the probability density function of regional (and more local) DSL values and 10-year return extreme values over the period 2000-2100 and compare these results to those obtained from a similarly forced non-eddying version of POP.  (Maltrud et al., 2008). The high spatial resolution captures the processes leading to mesoscale ocean eddies and provides a more detailed representation of the western boundary currents.

75
Specific details about the high-resolution model setup, such the treatment of the bottom topography, sea ice and river runoff, are described in the Weijer et al. (2012). The high-resolution model was optimised for use on the Cartesius supercomputer in Amsterdam (www.surfsara.nl) and about 3 model years are simulated per 24h using about 1000 cores.
The POP simulation was initialised from a 75 year spin-up simulation (Maltrud et al., 2008) under 80 the CORE-I climatology dataset (Large and Yeager, 2004) as atmospheric forcing. This initial condition is indicated here as the year 1950. Under a freshwater flux which is diagnosed from the last five years of this spin-up, the model displays only a very small drift over a 200 year control simulation (Bars et al., 2016). Here, the model was forced with monthly mean atmospheric forcing fluxes over the period 1950-2100, that were derived from simulations with the ECHAM5-OM1 model within the 85 ESSENCE (Sterl et al., 2008) project (see www.knmi.nl/∼sterl/Essence/). The used forcing fields are 10 m wind speed, downward flux of short wave and longwave radiation, 2 m temperature, humidity, precipitation, runoff, and the surface wind stress field. The atmospheric forcing fields are given on a global 1 • × 1 • grid and are interpolated to the curvilinear POP model grid. The outgoing heat and freshwater fluxes are computed within the model using bulk formulae. There is an initial adjustment  (Sterl et al., 2008)). We choose ensemble member #021 for which the high-resolution simulation is denoted by R 021 . In addition to this high-resolution simulation, a similarly forced simulation is performed with the low resolution POP version, indicated in the following by R low 021 . This non-eddying version has an average 1.0 • horizontal 100 resolution and 40 vertical levels (Weijer et al., 2012). The Gent and McWilliams (1990) scheme is used to represent eddy-driven tracer transports. Such a scheme is not needed in the strongly eddying version as these tracer transports are explicitly resolved.
The POP model directly computes the DSL, which can be decomposed into a mass redistribution term and a steric contribution. Because the freshwater flux is included into the model as a virtual 105 salt flux and the global mean of precipitation, evaporation and river runoff is zero, no mass induced global mean sea level changes can be represented. Due to the applied Boussinesq approximation, global mean steric sea level variations are not accounted explicitly during this study but this spatially independent contribution was computed from the model output (Greatbatch, 1994).
To demonstrate the performance of both versions of the POP model, we compare the DSL over the agrees well with observations, both for the mean (compare Fig. 1a and Fig. 1c) and the standard deviation (compare Fig. 1b and Fig. 1d). This shows that the model adequately determines the mean 115 ocean circulation, including the western boundary currents and also represents the eddy-induced variability. Differences with respect to observations appear due to the general overestimation of the modeled variability during this time period, which may be due to the prescribed low resolution of the atmospheric forcing and the lack of feedback of the atmosphere on the ocean variability. Differences in variability may also occur due to the higher horizontal resolution of the ocean model (0.1 • ) than 120 the altimetry dataset used (0.25 • ) as more small-scale features can be resolved. Regional differences in variability in the South Atlantic are caused by a too regular Agulhas ring formation rate in POP compared to observations (Bars et al., 2014).
In contrast, results for the low resolution simulation R low 021 shown in Fig. 1e  shows too low variability, in particular in the regions of the western boundary currents similar to many other non-eddying ocean model results (Bordbar et al., 2015). This weak variability also has consequences for the mean flow through the lack of representation of rectification processes causing, Monthly mean data are used for the analysis of changes in the mean and standard deviation (section 3.1), while daily data are used in the extreme value analyses (section 3.2).

135
In the POP simulation R 021 , global mean steric height increases by about 2.2mm/yr from year 2000 to 2100. As this signal is homogeneous over the Earth, it is not considered in the results below.
Largest changes in mean DSL between the periods 2081-2100 and 2001-2020 occur in the North Atlantic (Fig. 2a), in particular near the western part of this basin (Fig. 2c). There is a mean DSL decrease in the Atlantic and Pacific parts of the Southern Ocean, while mean DSL increases in the 140 Indian part of the Southern Ocean. The mean DSL increases in the eastern part of the North Atlantic basin and decreases in the center of the subpolar gyre. Large changes in DSL variability occur in the Agulhas retroflection region and near Drake Passage (Fig. 2b). The DSL variability decreases in the western North Atlantic, in the center of the subpolar gyre and slightly along the western boundary of the North Atlantic while it substantially increases in the eastern Atlantic (Fig. 2d).

145
In the POP simulation R low 021 , global mean steric height varies only by a few cm over the period 2000 to 2100 and again is not considered further. Regarding mean DSL patterns and amplitudes, the results of the low resolution simulation (R low 021 ), as shown in Fig. 2e, agree well with many other model studies using non-eddying ocean models (Landerer et al., 2007;Yin et al., 2009Yin et al., , 2010Bordbar et al., 2015). At first sight, the results also look similar to those for the R 021 simulation (compare 150 Fig. 2a and Fig. 2e). However, when regional details are considered, the results are different. The Southern Ocean basin contrast (Indian versus Atlantic/Pacific) is much stronger in the R 021 results.
The DSL change in the Northern Atlantic is more dipolar in the North Atlantic in the R low 021 results, with a large area south of Greenland with decreasing mean DSL. The change in DSL variability is, as expected, different in both models (compare Fig. 2b and Fig. 2f), in particular in western boundary 155 current regions. In the North Atlantic, (compare Fig. 2d and Fig. 2h), the changes in variability are less coherent in the Gulf Stream region and have larger amplitudes in the eastern part of the basin.
To explain the changes in DSL in the North Atlantic (Fig. 2c-d), the behavior of the Atlantic Meridional Overturning Circulation (AMOC) is shown in Fig. 3. The maximum AMOC at 26 • N decreases from about 20 Sv to about 5 Sv (red curve in Fig. 3a). The spatial pattern of the AMOC 160 is not changing but the North Atlantic Deep Water is shallowing by about 1000 m (Fig. 3b-c). The maximum strength of the AMOC at 35 • S decreases (blue curve in Fig. 3a) by more 60%. The decline in the AMOC causes a rise in mean DSL of up to 0.4 m near the North American continent, mostly because of a redistribution of ocean mass towards these regions (cf. Fig. 1a). The reduction of the AMOC also causes a northward shift of the latitude separation of the Gulf Stream. This result has also been found in the non-eddying model studies (Landerer et al., 2007) and previous strongly eddying model studies (Brunnabend et al., 2014). In addition, eastward shifts of the path of the Gulf Stream and North Atlantic Current occur. This is shown more clearly by the In the R 021 simulation, the global mean sea surface temperature (SST) rises by about 2 • C over the period 2000-2100 (Fig. 5a). Almost all ocean regions experience a warming and near the east coast of North America, there is a warming of up to 4 • C. However, in the subpolar gyre region of 180 the North Atlantic, SST is decreasing by more than 3 • C. SST remains almost unchanged over large regions in the Southern Ocean. This can be explained by the atmospheric forcing fields associated with the SRES-A1B scenario as they lead to changes in the radiative forcing between atmosphere and ocean. In addition, the decrease of the AMOC strength reduces the heat transport to the northern polar regions, which cools the upper ocean in the subpolar region. This leads to thermal contraction, 185 and a negative DSL change (Fig. 2a).
The DSL change in the Southern Ocean between the periods 2081-2100 and 2001-2020 ( Fig. 2a) is caused by a southward shift of the westerly winds. In addition, the westerly wind stress strengthens by about 0.03 Pa (Fig. 5b). The increase in zonal momentum flux accelerates the Antarctic Circumpolar Current and increases the northward Ekman transport that changes the slope of the isopycnal 190 surfaces in the South Atlantic (Yin et al., 2010). These effects cause changes in the water mass properties leading to steric contraction in the Southern Ocean and steric expansion in the region of the Agulhas return current (Yin et al., 2010), explaining the results in Fig. 2a.
The reduction of the AMOC also decreases the ocean-atmosphere temperature difference in the subpolar Atlantic region and hence leads to a reduction in the net ocean-atmosphere surface heat 195 flux, i.e. a reduced heat loss to the atmosphere (Fig. 5c). The cooling leads to reduced evaporation resulting in a freshening of the upper ocean (Fig. 5d), which causes a further slowdown of the AMOC as the AMOC is particularly sensitive to freshwater anomalies in this region (Smith and Gregory, 2009;Weijer et al., 2012). The changes in surface fluxes for the simulation R low 021 (not shown) are very similar as they are derived from the same atmospheric forcing fields, and are only slightly differently 200 affected by the ocean fields, compared to the R 021 simulation. Because the mechanism of deep water 6 Ocean Sci. Discuss., doi:10.5194/os-2016Discuss., doi:10.5194/os- -57, 2016 Manuscript under review for journal Ocean Sci. Published: 2 August 2016 c Author(s) 2016. CC-BY 3.0 License. formation is very different in the low-resolution model, the AMOC responds more mildly to changes in surface forcing than that in the high-resolution model (Weijer et al., 2012).

Regional PDF and Extremes
To determine an estimate of the Probability Density Function (PDF) of DSL we show histograms of The changes in each PDF for the R 021 simulation for the different regions and locations are plotted in Fig. 7 with the blue histogram being the future PDF. The variance of DSL decreases in 215 mid-Atlantic region 1 (cf. Fig. 2b,d), which is seen by the shift of the PDF to the left (Fig. 7a). This also leads to a reduction of the highest DSL extremes by more than 10 cm. In region 2 (western North Atlantic) DSL is mainly driven by mean changes due to steric effects and the mass redistribution and hence the PDF shifts to the right (Fig. 7b). In the eastern North Atlantic (region 3), the variance of the DSL increases (cf. Fig. 2b,d) due to the changes in the pathways of eddies causing the changes 220 in EKE (Fig. 4c,d). This leads to a rightward shift of the PDF by about 10 cm in this region (Fig.   7c).
Changes in the pathways of eddies are also important when considering local DSL extremes. The Azores are located in a region of slightly decreased variability (Fig. 1b,d) due to reduced eddykinetic energy in this region, shifting the PDF slightly to the left (Fig. 7d). Near the Bermuda Islands 225 the shift in the ocean currents leads to lower probabilities of higher sea level extremes. (Fig. 7e).
The most interesting result, however, is shown in Fig. 7f for the coast near Lisbon. Due to the shift in the Gulf Stream and North Atlantic Current one would expect increased probabilities for high DSL values in this region. However, because these currents are not only shifted but also reduced in strength almost no changes in DSL extremes can be identified (Fig. 7f).

230
The changes in the PDFs for the R low 021 simulation show quite a different behavior than those in the R 021 simulation for most regions and locations. While the relative shift in the mean is comparable for both models in the regions 1 and 2 (Fig. 8a,b), the amplitude is much smaller for R low 021 . For region 3 (Fig. 8c), the PDF has bimodal characteristics and hardly changes under climate change, in contrast to the change in the R 021 simulation (Fig. 7c). The PDF change for the Azores is the influence on ocean eddy paths in the R low 021 simulation (Fig. 4b). The PDF of the other two locations (Fig. 8e,f) show the same behavior as in the R 021 simulation.
The extreme DSL values for a specific return time are computed here in the same way as in Brunnabend et al. (2014). i.e. from a fit of parameters in a Generalized Extreme Value (GEV) dis-240 tribution (Coles, 2001). The extreme DSL values at a return time of 120 months (10 years) over the period 2001-2020 and their changes over the different 20 year periods (2081-2100 and 2001-2020) of the R 021 simulation are shown in Fig. 9. Over the period 2000-2020 higher extreme sea levels occur in regions of high variability, i.e. in regions of the major current systems such as the Gulf Stream and the Agulhas Current (Fig. 9a,c). Therefore, the regional pattern of changes in extreme sea lev- The same extreme sea level values are shown in Fig. 10 for the R low 021 simulation. The amplitude of these extremes is much smaller, in particular in western boundary current regions (Fig. 10a) and in the Gulf Stream region (Fig. 10c). The low-resolution ocean model simulation leads to different 255 extreme sea level projections in the northern North Atlantic (in particular, in the Labrador Sea and Barents Sea) than for the R 021 simulation. The sign of the change in sea level extremes is also different in the Caribbean Sea. This shows the importance of including an explicit representation of eddy processes into an ocean model when looking at regional projections of DSL.

260
In this paper, we considered a single realisation of future dynamic sea level (DSL) changes using a strongly eddying ocean model forced by atmospheric fields according to an SRES A1B scenario.
The results show that changes in local and regional PDFs (between the periods 2001-2020 and 2081-2100) are mainly due to changes in DSL variability on short time scales and therefore related to changes in the ocean eddy field. This can be deduced from both the changes in the eddy kinetic 265 energy of the ocean surface velocity field and from a comparison of DSL changes in a non-eddying version of the same model. In the high-resolution model simulation, the changes in eddy pathways are caused by a strong decrease of the AMOC with simultaneous eastward shifts in the path of the Gulf Stream and North Atlantic Current.
Our main result is that the pattern of 10-year return time DSL extremes (as shown in Fig. 9) are 270 determined by changes in the ocean eddy field (Suzuki et al., 2005;Brunnabend et al., 2014). In 8 Ocean Sci. Discuss., doi:10.5194/os-2016Discuss., doi:10.5194/os- -57, 2016 Manuscript under review for journal Ocean Sci. Published: 2 August 2016 c Author(s) 2016. CC-BY 3.0 License. some regions of the globe these extreme DSL values can be up to 0.5 m which of the same order of magnitude as the mean DSL change. This shows the importance of internal ocean variability for regional extreme sea levels, not only on the longer time scales (Bordbar et al., 2015), but also on the shorter time scales (Firing and Merrifield, 2004).

275
Low resolution ocean/climate models are not capable of accurately representing these changes in extreme sea levels. Some low resolution model studies do capture a shift in ocean currents in case of a declining AMOC (Landerer et al., 2007;Yin et al., 2010Yin et al., , 2009Pardaens et al., 2011;Kienert and Rahmstorf, 2012). However, the model resolution is not resolving DSL variability caused by ocean eddies, as the parameterisation of eddies in these models only affects the heat and salt transport in the 280 models. Already using an eddy-permitting ocean/climate model (with a 0.25 • horizontal resolution) indicated the importance of resolving ocean eddies to accurately estimate future sea level variability (Suzuki et al., 2005), but in these models the western boundary currents usually do not have a correct separation behaviour.
There are several caveats in this model study which may modify the results quantitatively but 285 which do not affect the main message of this paper that strongly eddying models are important for regional future sea level change projections. First, the AMOC in the R 021 POP model simulation appears to be very sensitive to freshwater anomalies and hence the scenario here may be quite an extreme one. Second, the usage of an ocean-only model with mixed boundary conditions below sea-ice regions and atmospheric forcing fields from a climate model is restricting the capabilities of 290 the model in simulating the coupled ocean-atmosphere interactions occurring in reality. However, it is expected that shifts in the ocean eddy fields would also occur in coupled models with strongly eddying ocean model components. Third, the model does not simulate many other processes causing regional and coastal sea levels changes (e.g. GIA, and gravity). Many of these processes would only effect the mean DSL values and not its variability. Hence, as a first approximation, these sea level 295 changes can be added to the mean DSL values (Slangen et al., 2012(Slangen et al., , 2014 determined here. Finally, we studied only one realization of the model here and it would be better to use an ensemble of simulations (Bordbar et al., 2015) to determine the effect of ocean initial conditions and to have better statistics on the extreme DSL values. The latter is still hardly feasible with the current computational capabilities. 300 We conclude from the results that when developing plans for adapting to future changing sea level, not only mean regional changes should be considered, although they may be substantial. Also the changes in variability should be accounted for, as with higher variability the probability of sea level extremes may increase. This in particular holds for the North Atlantic region where many areas are vulnerable to sea level rise.