Ocean model resolution dependence of Caribbean sea-level projections

Sea-level rise poses severe threats to coastal and low-lying regions around the world, by exacerbating coastal erosion and flooding. Adequate sea-level projections over the next decades are important for both decision making and for the development of successful adaptation strategies in these coastal and low-lying regions to climate change. Ocean components of climate models used in the most recent sea-level projections do not explicitly resolve ocean mesoscale processes. Only a few effects of these mesoscale processes are represented in these models, which leads to errors in the simulated properties of the ocean circulation that affect sea-level projections. Using the Caribbean Sea as an example region, we demonstrate a strong dependence of future sea-level change on ocean model resolution in simulations with a global climate model. The results indicate that, at least for the Caribbean Sea, adequate regional projections of sea-level change can only be obtained with ocean models which capture mesoscale processes.

Scientific RepoRtS | (2020) 10:14599 | https://doi.org/10.1038/s41598-020-71563-0 www.nature.com/scientificreports/ global climate model is resolution dependent 15,16 . From globally eddy-resolving ocean model simulations 17,18 , it is known that explicitly capturing eddies leads to substantially different DSL responses compared to lowerresolution models, for example through the modification of boundary currents. Regional sea-level projections based on global climate models which parameterise mesoscale processes (such as in CMIP6) therefore miss relevant physics 19 affecting the projected sea-level rise at the end of the century.
To study the effect of ocean model resolution on DSL and SDSL projections, we here analyse the Caribbean region. This is a typical region where large effects can be expected as it lies downstream of a strong western boundary current, the North Brazil Current (NBC). This current is characterised by the shedding of mesoscale ocean eddies which strongly affect the downstream region, the Caribbean Sea. The NBC is also part of the larger Atlantic Meridional Overturning Circulation (AMOC) which is expected to weaken under climate change 20,21 . Both the large-scale ocean circulation and the mesoscale eddies are expected to be much better represented in a high-resolution ocean model compared to a low-resolution one and here below we explore the consequences for regional (stero)dynamic sea-level projections.

climate model simulations
In order to systematically investigate the effect of ocean model resolution on simulated Caribbean sea-level rise, we performed simulations with the Community Earth System Model (CESM). The CESM is a fully-coupled climate model with a volume conservation constraint for the ocean component. The CESM has no dynamic ice sheet model and hence the effects of mass loss by glaciers and of the Greenland-and Antarctic ice sheets are not captured. The high-resolution version of CESM has an ocean component with a 10 km ( 0.1 • ) horizontal resolution capable of capturing the development and interaction of mesoscale ocean eddies 11 and an atmosphere component with a horizontal resolution of 50 km ( 0.5 • ). The ocean (atmosphere) component of the low-resolution version of the CESM has a horizontal resolution of 100 km (125 km). The ocean component of this low-resolution model does not capture mesoscale processes.
With the high-resolution version, we first performed a 200 years spin-up with seasonally varying yearlyrepeated forcing conditions (aerosols, solar insolation) of the year 2000 (with a pCO 2 level of 369 ppmv). After the 200 year spin-up period, we branched two simulations, one in which the spin-up was further extended (HR-CESM Control) and one in which the atmospheric pCO 2 increases by about 1% each year (HR-CESM) for 101 years (Supplementary Figure S1). The HR-CESM simulation is the first of its kind due to its extremely high computational costs. The low-resolution version of the CESM was spun-up with a present-day configuration (similar to the HR-CESM) for 1,200 years and we branched a control simulation (LR-CESM Control) and a simulation with about 1% pCO 2 increase each year (LR-CESM) for 101 years. Both the control simulations are almost free of any temperature trends over the upper 1,000 m of the ocean compared to the HR-CESM and LR-CESM simulations (Supplementary Figure S2). More details of the CESM simulations are provided in "Methods" section.
The majority of the CESM output is stored as monthly-averaged quantities. These monthly-averaged quantities are converted to yearly averages which are used in the analyses below. A limited number of quantities are stored as daily averages which are used in generalised extreme value analysis (see "Methods").

Sterodynamic sea-level trends
We first determined the trends in the yearly-mean SDSL field (indicated by η , see "Methods") over the 101-year period (2000-2100) for the HR-CESM and LR-CESM simulations. The local SDSL change consists of two components: ocean dynamic sea-level change (indicated by η M , DSL change) and global-mean thermosteric sealevel rise (indicated by η g S ) 22 . The η g S of the HR-CESM and LR-CESM are corrected for any drift in their control simulations and are fairly similar when comparing the HR-CESM and LR-CESM (see Table 1; Supplementary Figure S7). The largest contribution ( > 80% ) to the total η g S originates from the upper 1,000 m of the ocean. Both simulations show a positive (and significant) SDSL trend over the displayed Caribbean region (Fig. 1a,b) and differences become more pronounced when these trends are normalised by the trend in η g S (Fig. 1c,d). The fastest SDSL rise occurs near the southern and western boundary of the Caribbean Sea in both simulations. However, the normalised SDSL trends in these parts of the Caribbean Sea are above average and below average with respect to the η g S trend for the HR-CESM and LR-CESM, respectively. The differences in the (normalised) SDSL trends between the HR-CESM and LR-CESM are related to DSL changes. The local DSL ( η M ) trend patterns are shown in Supplementary Figure S3a,b for the HR-CESM and LR-CESM, respectively, and there are large differences in both amplitude and sign (see also Supplementary Figure S4a). We also determined the local steric sea-level (indicated by η S ) change and the trends are shown in Supplementary Figure S3c,d. The local η S is corrected for any drift using the control simulations. The local η S trend is positive in the displayed region and the slowest η S trends are found on continental shelf for both CESM simulations. The η S trends in the Caribbean Sea and surroundings are faster in the HR-CESM compared to the LR-CESM (Supplementary Figure S4b).
The differences in the sign (above or below average) of the normalised SDSL trends across the Caribbean Sea (north-south or east-west) is also present in altimetry observations (Supplementary Figure S5). This dipole pattern in the normalised SDSL trends is related to the NBC and Caribbean Current, which typically determine the DSL distribution. The oceanic fronts and currents are much better represented in HR-CESM than LR-CESM 23 , leading to a different response in DSL in both simulations, of which the HR-CESM is more in agreement with observations.
Note that the observation record is much shorter (25 years) compared to the CESM simulations (101 years). Therefore we determined the normalised SDSL trend over a shorter period for the CESM simulations (Supplementary Figure S6). Here we present results of the normalised SDSL trends over both regions 1 and 2 ( Fig. 1), which are situated west (region 1) and east (region 2) of the western boundary current. After model year 2050, Scientific RepoRtS | (2020) 10:14599 | https://doi.org/10.1038/s41598-020-71563-0 www.nature.com/scientificreports/ the normalised SDSL trends have similar magnitudes and signs as the normalised SDSL trends over the entire period (2000-2100) and this results is robust for varying model initial year and the period over which the trends are determined. Before model year 2050, the magnitude and sign of the normalised SDSL trends are sensitive to the chosen period, which is likely related to natural variability as the simulations are initiated from a control simulation. The normalised SDSL centennial trends are representative for the shorter 25-year period (same length as observations) normalised SDSL trends in the second half of the CESM simulations, when the simulations are not in equilibrium. Both the present-day climate and the global-mean sea level are not in equilibrium either 6,8 . It is possible that observations shows persistent above-average and below-average sea-level trends for region 1 and region 2, respectively, due to natural variability, but the observational record is too short to falsify this hypothesis. The time evolution of the DSL ( η M ) and SDSL ( η ) averaged over both regions 1 and 2 are shown, together with the global-mean thermosteric sea-level rise ( η

changes in the large-scale ocean circulation
Long-term variations in the η M fields are related to variability in the large-scale ocean circulation. For example, a weaker western boundary current results in a decrease in the zonal pressure gradient according to the geostrophic balance. Locally, the decreased pressure gradient leads to an increase (decrease) of η M west (east) of the western boundary current.
One way to determine changes in the large-scale circulation (which affects the η M fields) is by analysing the vertically integrated flow, represented by the barotropic streamfunction (BSF). The BSF is computed by meridionally integrating the zonal component of the vertically integrated flow, see "Methods". The time-mean BSF fields at the beginning (model years 2000-2029) of both simulations are shown in Fig. 2a,b. The large-scale patterns of the subtropical gyre and subpolar gyre in the North Atlantic are represented in both CESM simulations. However, there are some notable differences in the BSF such as the Gulf Stream, which is better represented in the HR-CESM 24 as reflected by the sharper BSF gradient (around 40 • N , see Fig. 2a) compared to the LR-CESM (Fig. 2b).
The difference between the BSF over the years 2070-2099 and the years 2000-2029 is shown in Fig. 2c,d. Both the subtropical gyre and the subpolar gyre become weaker over time, with the largest negative BSF anomalies and meridional velocity anomalies occurring near western boundary currents (see also Supplementary Figure S8). Note the relatively large BSF anomalies in HR-CESM near 40 • N . Small changes in the path of the Gulf Stream 25 result in relatively large changes in BSF due to the sharp BSF gradients in that region. In LR-CESM, the BSF decreases almost over the entire Caribbean Sea and the BSF anomaly has a similar pattern as the normalised SDSL trends (Fig. 1d). We find significant correlations between the yearly-averaged η M and BSF fields in the Caribbean Sea over the 101-year period in the HR-CESM Control and LR-CESM Control (Fig. S9). www.nature.com/scientificreports/ Long-term changes in BSF affect the η M fields according to geostrophic balance. The large-scale reduction in the BSF over time could in principle be related to changes in the vorticity input by the wind (i.e., the wind stress curl), but we find no substantial differences in the wind stress curl between the early and later period of the simulations (compare contours in Fig. 2a,c and Fig. 2b,d). We also did not find any reduction of the wind stress curl at the north coast of South America which potentially could explain the positive η M trends in this part of the Caribbean Sea 26 .
A change in the Atlantic Meridional Overturning Circulation (AMOC) 20,21,27 can also affect the North Brazil Current and Caribbean Current strength 28,29 , since both currents are part of the northward branch of the AMOC. We find a significant decrease in the AMOC strength (see "Methods"), as seen through the insets in Fig. 2a,b and Supplementary Figure S10. The southward branch of the AMOC [(i.e. the North Atlantic Deep Water (NADW), see "Methods"] also decreases over time (insets in Fig. 2c,d, Supplementary Fig. S10). Hence changes in Atlantic Ocean circulation captured by the BSF, in particular those in the NBC, are mainly related to a reduction of the AMOC strength and consequently affect the η M fields in the Caribbean Sea.

future sea-level extremes
A prominent feature of the NBC is its retroflection which sheds off multiple anti-cyclonic eddies per year. These NBC eddies propagate along the background flow towards the Lesser Antilles and sometimes (partially) enter the Caribbean Sea 30 . To determine the effect of the NBC eddies (which are shed by the NBC retroflection) on the DSL, we analyse the local maxima of η M . The η M signature of the NBC eddies is partly filtered out in the monthly-averaged (or yearly-averaged) η M fields due to the time averaging. Therefore, we analyse daily-averaged η M fields over the 101-year period for both HR-CESM and LR-CESM simulations. Using the daily-averaged η M fields, we determine the local monthly maximum of η M (indicated by η max M ) for both simulations (insets in Fig. 3a,b). The cross here indicates the location of the maximum value, η Max M , over the NBC outflow region year event is only 1.5 cm ( −6% ). We found a similar decrease in η Max M when taking the yearly maxima (removing the seasonal cycle) but the yearly maxima GEV fits are less robust due to fewer data points (30 instead of 360).
Similar to the monthly η max M fields, we also determined the yearly η M maximum fields (retained from daily averages) and trends over the 101-year period for the simulations (see Supplementary Figure S12a,b). In the NBC outflow region, this trend in the yearly η M maximum is three to four times larger than the trend in the yearly-averaged η M (cf. Supplementary Figure S3c,d) for the HR-CESM (Supplementary Figure S12c). For the LR-CESM, the trend in the yearly η M maximum is only twice the magnitude of the yearly-averaged η M trend in the NBC outflow region (Supplementary Figure S12d).
Apart from the η M changes, we also find a weakening of the horizontal surface velocities in the NBC transport region (Supplementary Figure S13). Both HR-CESM and LR-CESM simulations show a path of significant ( 95%-confidence level) negative velocity trends from the NBC retroflection to the Lesser Antilles. This path of decreasing velocities in the LR-CESM is much broader compared to that in the HR-CESM.
To summarise, the differences in η M extremes and horizontal velocities in the NBC outflow region between the HR-CESM and LR-CESM simulations are related to the horizontal resolution of the models. In the LR-CESM simulation, the NBC retroflection and shedding of NBC eddies are parameterised and the relevant mesoscale dynamics in this region is not captured 11 . At least in the NBC outflow region, DSL projections based on the LR-CESM can result in overestimations in the 1:5 year moderately extreme event (over a period of 100 years) compared to the HR-CESM, with about a 6 cm (factor five) difference.

Sterodynamic sea-level trends in CMIP6 models
Most models participating in the Coupled Model Intercomparison Project phase six (CMIP6, see "Methods") do not explicitly capture mesoscale processes as they have a similar horizontal resolution as in the LR-CESM simulation (i.e., 1 • ). An analysis of the CMIP6 output (under the 1% pCO 2 increase scenario) shows an increase in η g S and a decrease in the AMOC strength over the simulated period, which is a similar response as in our CESM simulations (Supplementary Figure S14). In the CMIP6 models, there is a wide range in SDSL trends in the Caribbean Sea over a 101-year period (Supplementary Figure S15). The horizontal velocity trends (Supplementary Figure S16) are similar to our LR-CESM simulation, even for the eddy-permitting models (i.e. 0.25 • ). Figure 4b indicates that, from comparing the normalised SDSL trends (with respect to the η g S trend) over regions 1 (south-west Caribbean Sea) and 2 (north-east of Caribbean Sea), the HR-CESM has the same normalised SDSL signs as observations for both regions (i.e. in the same quadrant). Note that the observed normalised sea-level trends (by 3 mm year −18 ) include all contributions of sea-level rise, so it is not the 'pure' SDSL trend. The contribution of melt by the Greenland-and the Antarctic ice sheets is distributed homogeneously over the Caribbean region 9 and will not affect the sign of the observed normalised sea-level trends (it will change the numerical values). Ten (Two) out of the fifteen CMIP6 models show positive (negative) DSL trends for both regions (Fig. 4a). Only three CMIP6 models have the same normalised SDSL signs as observations for both regions (Fig. 4b). The normalised SDSL CMIP6 mean and standard deviation is 1.27 ± 0.22 and 1.05 ± 0.22 for region 1 and region 2, respectively. The normalised sea-level trend in observations is 1.14 and 0.82 for region 1 and region 2, respectively. Note that the observed trends are determined over a shorter period compared to the centennial CESM and CMIP6 trends, as already discussed in Supplementary Figure S6. The CMIP6 mean does not have the same SDSL sign as observations for region 2. The normalised SDSL of the HR-CESM (as well as the LR-CESM) lies outside the CMIP6 standard deviation for region 2. This mismatch for region 2 is related to the misrepresentation of the effects of NBC eddies, which are not captured by the CMIP6 models. From these preliminary CMIP6 model results, the largest biases in SDSL are found in region 2. Heat uptake by the (upper) ocean causes global-mean thermosteric sea-level rise ( η g S ), which is the main contribution of local sterodynamic sea-level rise (at least 75% ). Although the effects of ice-sheet mass loss and isostatic adjustment are not taken into account in any of the presented model simulations, these are expected to be fairly homogeneous on the regional scale 9,32 and can be incorporated in the global-mean thermosteric sealevel rise (e.g. in η g S in this study). The changes in η M are mainly related to a reduced Atlantic Meridional Overturning Circulation (AMOC) and the associated changes in the North Brazil Current (NBC). In the HR-CESM simulation, the reduced strength of the NBC results in a decreased eddy activity and strength, resulting in weaker η M extremes. For the 1:5 year moderately extreme event in the maximum of η M over the NBC region, there can be a factor 5 difference in DSL values between HR-CESM and LR-CESM.

conclusions and implications
In the Caribbean Sea and surroundings, mesoscale processes (such as the NBC retroflection and eddies) are important for regional sea-level projections. Adequate Caribbean sea-level projections can only be obtained in ocean models which capture mesoscale oceanic processes. Consequently, most CMIP6 models as well as the LR-CESM are not fit for purpose in making such projections. The higher resolution ocean models (AWI-CM-1-1-MR  Table 1) for the HR-CESM, LR-CESM, the CMIP6 output and (in b)] AVISO for the two different regions 1 and 2 (cf. Fig. 1). In (b), the AVISO determined observations (the total sea level, η ) are normalised with a value of 3 mm year −1 . The colour-coded boxes in the quadrants indicate above-average (red) and below-average (blue) SDSL trends with respect to the η g S trend for the two regions.
Scientific RepoRtS | (2020) 10:14599 | https://doi.org/10.1038/s41598-020-71563-0 www.nature.com/scientificreports/ and CNRM-CM6-1-HR) of the CMIP6 protocol with a 0.25 • ocean model resolution 33 , are not noticeably better compared to the coarse CMIP6 models in the Caribbean region. We have shown that the effect of DSL can locally have a large contribution to the sterodynamic sea-level change. The DSL is affected by mesoscale processes everywhere, but particularly in eddy-active regions such as near western boundary currents. Hence, although we only showed results for the Caribbean, the effects of ocean model horizontal resolution will be important for sea-level projections in other regions and needs to be further investigated. Due to the high-computational costs, only one realisation of the high-resolution model was available, yet the results here indicate that an ocean model with a 0.1 • horizontal resolution agrees much better to observations and this appears to be the minimum resolution required to capture an adequate dynamic sea-level field.
For the Caribbean region, the good news is that changes in dynamic sea-level extremes in the high-resolution model projects a much lower sea-level rise than the low-resolution model, in particular for the many island areas in the eastern part of the region. High-resolution climate models are crucial for a better understanding of future sea-level rise and sea-level extremes, and for planning future investments to adapt to effects of climate change in the Caribbean region and elsewhere over the globe. Supplementary Figure S1 shows some key variables for the spin-up period (similar results as in 36 ), the control simulation and the forced simulation for the high-resolution model. For the control simulation, the global mean (2 m) surface temperature is almost constant and the radiative imbalance is slightly positive. The Gregory plot (Supplementary Figure S1c) also shows the equilibration of the control simulation and the deviation from the equilibrium for the forced simulation. The upper 700 m global ocean heat content (Supplementary Figure S1d) is still adjusting, but relatively small trends (compared to the spin-up) occur over the last 100 years. The deep ocean fields take a much longer time to equilibrate. In the forced simulation, the surface temperature, radiative imbalance and ocean heat content start to deviate from the control simulation due to the increase of atmospheric CO 2 concentration.

Methods
In addition to the high-resolution CESM simulations, a companion simulation was conducted at a lower resolution (using CESM version 1. CMIP6 data. We use results from the latest release of the Coupled Model Intercomparison Project phase six (CMIP6) and compare these to the output of our CESM simulations (Table 1). We analysed the preliminary model output of the CMIP6 in which atmospheric CO 2 levels increase each year by 1% . Note that each model of the CMIP6 is initiated from the pre-industrial (year 1,850) control simulation. Only models which conserve ocean volume (as the CESM) are considered here (variable name 'zos'). We analyse the monthly averaged η M , horizontal velocity, temperature and salinity fields of the first 101 model years, as is done for the CESM output. The model output of the AWI-CM-1-1-MR is provided on an unstructured grid. Before analysing the output of this CMIP6 model, the data is interpolated onto the HR-CESM grid.
Significance of the linear trends. The trends are derived from a linear least-square fit to the yearly-averaged time series. The significance of each trend is determined following the procedure outlined in 37  The Student-t value is the ratio between the linearly fitted trend and the standard error. Using the reduced degrees of freedom and the two-sided critical Student-t values, one can determine the significance of having a trend different from zero (the null hypothesis).
Sterodynamic sea level. The local sterodynamic sea level (SDSL) consists of two components 22 . The first component is the dynamic sea level ( η M ) and is the height of the ocean surface above the geoid. The η M fields are part of the standard output in the CESM simulations (variable 'SSH') and CMIP6 models (variable 'zos'). The η M fields have a global mean of about 0, if the global mean was non zero we subtracted uniformly the global mean from the η M fields. The next component is the global-mean thermosteric sea-level rise ( η g S ), which is determined from postprocessing the model output 38 . First, we determined the local steric sea-level ( η S ). The contribution of both thermal changes and haline changes is determined as the full-depth integral over the specific volume anomaly 39 , The temperature, salinity and pressure dependency are taken into account while determining the density, and ρ 0 = 1, 028 kg m −3 . The steric sea-level change is expressed as an anomaly with respect to the initial value of the first model year. The global-mean thermosteric sea-level ( η g S ) rise is determined by globally averaging η S . This procedure is done for both the CESM simulations as the CMIP6 models. One can use the variable 'zostoga' instead of η g S for the CMIP6 models, but this is not used in this study for comparison with the CESM simulations. Eventually, hence the local sterodynamic sea level ( η ) is determined from The trend over the 101-year period for η M and η S are shown in Supplementary Figure S3. We corrected for any drift in η S and η g S using the control simulations.
Barotropic streamfunction. The barotropic flow is defined as the full depth integral of the horizontal velocity: (1) Table 1. Overview of the CMIP6 and CESM models with the resolution of the ocean grid. The global-mean thermosteric sea-level ( η  Generalised extreme value. The η M extremes in the NBC outflow region are analysed using GEV (Generalised Extreme Value) analysis. First, we retained the monthly η M maxima ( η Max M ) from daily averaged η M . This is the block maxima GEV approach where each month contains the η M maximum. The distribution of η Max M may vary in time. Therefore different periods of n months' length are retained from the full time series, under the assumption that these periods in time are stationary and stationary GEV analysis can be conducted 41 . Using all the data points over a particular period, a GEV fit is made to the data points using the following expression as the cumulative distribution function: where µ , σ , ξ are the location, scale and shape parameter, respectively. A special case of the distribution is where ξ → 0 , the GEV distribution approaches a Gumbel distribution. The probability of occurrence (p) using the distribution of G n (z) is related to return time ( τ ) of such an event: The value of η Max M related to this probability ( z p ) can be obtained by: Note that detrending the periods of time (linear or quadratic) might result in biases. Therefore, the data on which GEV is applied is directly fitted to the G n (z) distribution. We have chosen to fit the G n (z) distribution using periods of 360 months' length (30 years). Shorter sections (120 months) provided similar results, but might have insufficient data to make a reasonable fit.

code availability
All the model output is analysed with Python 2.7.9. The maps are generated using the Basemap package in Python. All figures are prepared with Python 2.7.9. Part of the Python code as well as processed CESM and CMIP6 model output can be accessed at https ://githu b.com/Renev anWes ten/SR-Carib bean.