Yin et al. (2020)

Title:: Response of Storm-Related Extreme Sea Level along the U.S. Atlantic Coast to Combined Weather and Climate Forcing
Keywords:: Meridional overturning circulation, Extreme events, Sea level, Storm surges, Climate change, Climate models
Corresponding author:: Jianjun Yin
Citation:: Yin, J., Griffies, S. M., Winton, M., Zhao, M., & Zanna, L. (2020). Response of Storm-Related Extreme Sea Level along the U.S. Atlantic Coast to Combined Weather and Climate Forcing. Journal of Climate, 33(9), 3745–3769. doi: 10.1175/jcli-d-19-0551.1
URL:: https://journals.ametsoc.org/view/journals/clim/33/9/jcli-d-19-0551.1.xml

Abstract

Storm surge and coastal flooding caused by tropical cyclones (hurricanes) and extratropical cyclones (nor’easters) pose a threat to communities along the Atlantic coast of the United States. Climate change and sea level rise are altering the statistics of these extreme events in a rather complex fashion. Here we use a fully coupled global weather/climate modeling system (GFDL CM4) to study characteristics of extreme daily sea level (ESL) along the U.S. Atlantic coast and their response to global warming. We find that under natural weather processes, the Gulf of Mexico coast is most vulnerable to storm surge and related ESL. New Orleans is a striking hotspot with the highest surge efficiency in response to storm winds. Under a 1% per year atmospheric CO2 increase on centennial time scales, the anthropogenic signal in ESL is robust along the U.S. East Coast. It can emerge from the background variability as soon as in 20 years, or even before global sea level rise is taken into account. The regional dynamic sea level rise induced by the weakening of the Atlantic meridional overturning circulation facilitates this early emergence, especially during wintertime coastal flooding associated with nor’easters. Along the Gulf Coast, ESL is sensitive to the modification of hurricane characteristics under the CO2 forcing.

Introduction

The U.S. Atlantic coast (including both the East Coast and the Gulf of Mexico coast) is an active region for tropical and extratropical storms. Geographically, this densely populated coastal region is surrounded by a broad (100–300 km) and shallow (<100 m) continental shelf (Fig. 1a), making it particularly vulnerable to severe storm surge and associated socioeconomic damages and life loss. Notable examples include Hurricane Katrina in 2005 (Fritz et al. 2007) and Superstorm Sandy in 2012 (Sobel 2014), as well as more recent wintertime coastal flooding caused by “bombogenesis” (Buell 2018).

Figure 1: Geometry and bathymetry along the U.S. Atlantic coast. (a) The U.S. Atlantic coastline in nature is highlighted by the red color. Large coastal cities are marked from north to south and west: Halifax, Boston, New York, Baltimore, Charleston, Miami, Tampa, New Orleans, and Houston/Galveston. (b) Representation of coastal geometry and bathymetry (km) in CM4. Coastal ocean grid boxes in red, green, and purple indicate the NE (north of Cape Hatteras), SE (Cape Hatteras to the south tip of Florida), and GOM regions. The color scale uses 100-m intervals for 0–500-m depth and 500-m intervals below 500 m.

The twenty-first-century outlook of storm surge often invokes the “noise + trend” model, namely, global sea level rise (SLR) will lead to elevated storm surge (USGCRP 2017). However, the real world situation is more complex. Storm surge critically depends on such storm characteristics as intensity, frequency, size, path, translational speed, and landfall angle (Simpson 1974; Weisberg and Zheng 2006; Irish et al. 2008; Rego and Li 2009; Hall and Sobel 2013). In addition to chaotic/stochastic “noise” processes, the generation, development, and propagation of tropical cyclones (TCs) and extratropical cyclones (ECs) over the North Atlantic and North America are influenced by El Niño–Southern Oscillation (ENSO) (Hirsch et al. 2001; Donnelly and Woodruff 2007), the North Atlantic Oscillation (Elsner et al. 2000; Ezer and Atkinson 2014), the Atlantic multidecadal variability (Zhang and Delworth 2006), anthropogenic greenhouse gas and aerosol forcing (Mann and Emanuel 2006; Lin et al. 2012; Little et al. 2015; Garner et al. 2017; Rahmstorf 2017; Marsooli et al. 2019), and other factors. These factors mutually interact and modify sea surface temperature (SST), ocean heat distribution, vertical wind shear, meridional temperature gradient, large-scale oceanic and atmospheric circulation, and regional and global sea level. Given their distinct spatiotemporal scales and possible opposing effects on storms and storm surge, it remains scientifically challenging to study sea level extremes and their impact along the U.S. Atlantic coast in the face of natural and anthropogenic climate variability and change.

Here we address this challenge using a fully coupled global climate model (CM4) recently developed at the Geophysical Fluid Dynamics Laboratory (GFDL) of NOAA. Under the protocol of phase 6 of the Coupled Model Intercomparison Project (CMIP6) (Eyring et al. 2016), a series of climate change experiments have been performed with CM4. With these simulations, we focus on weather–climate interactions and their combined effect on storm-related extreme daily sea level (ESL) along the U.S. Atlantic coast. The paper is organized as follows. Section 2 describes the model and daily sea level analysis methods. Section 3 evaluates the model performance in sea level simulations. Section 4 presents characteristics and statistics of ESL under natural weather processes. Their response to CO2 forcing is described in section 5, followed by the conclusions and discussion of model limitations toward the end.

Model, data, and methods

The GFDL CM4 climate model

CM4 is the latest generation of the climate models developed and used at GFDL (Held et al. 2019). The atmospheric model (AM4.0) (Zhao et al. 2018a,b) adopts finite-volume cubed-sphere dynamical core with 96 (~1.0° grid spacing) or 192 (~0.5° grid spacing) grid boxes per cube face. It has 33 vertical levels and the model top is located at 1 hPa. The model incorporates updated physics such as a double-plume scheme for shallow and deep convection and single-moment cloud microphysics. Due to improvements in model resolution, physics, and dynamics, CM4 can better simulate strong TCs and ECs over the North Atlantic and North America with hurricane-force winds, reasonable storm tracks, seasonal cycle, and interannual variability (Zhao et al. 2018a), although the strongest (e.g., category 4 and 5) hurricanes are not simulated (Fig. 2).

Figure 2: Simulations of TC and EC in the long-term piControl of CM4. (a) Global map of TC tracks. (b) TCs over the North Atlantic. Black contours (hPa) indicate the mean subtropical high during June–October. (c) Global map of EC tracks. (d) ECs over the North Atlantic and North America. Black contours (°C) indicate the near-surface temperature and its gradient during December–February. The color scale denotes daily winds (m s−1) associated with storms.

We use a Lagrangian approach and the 6-h data for detecting and tracking TCs and ECs in the CM4 simulations (Fig. 2) (Zhao et al. 2009). Among multiple criteria, TCs should have a maximum surface wind speed of at least 14 m s−1 and their trajectories must last at least 3 days. In addition, TCs should have a warm core of at least 1°C above the surrounding temperatures between 300 and 500 hPa, thus allowing TCs to be distinguished from ECs and other storms. The relative vorticity at 850 hPa in TCs should be greater than 1.6 × 10^{−4} s^{−1}. We use the sea level pressure field to identify ECs that should have a maximum wind speed of 25 m s−1 and last for at least 3 days. Detailed algorithms and codes for detecting and tracking storms can be found at https://www.gfdl.noaa.gov/tstorms/.

The oceanic model of CM4 is based on the Modular Ocean Model version 6 (MOM6). It uses the arbitrary Lagrangian–Eulerian algorithm in the vertical to allow for the combination of different vertical coordinates including geopotential, isopycnal, and terrain following. The model adopts the C-grid stencil in the horizontal and is configured on a tripolar grid. It has a 0.25° eddy-permitting horizontal grid spacing (~20 km at midlatitudes) and 75 hybrid vertical layers down to the 6500-m maximum bottom depth. On the shelf, the vertical grid spacing can be as fine as 2 m. The ocean model configuration used here as part of CM4 is documented by Adcroft et al. (2019).

MOM6 roughly resolves important bays and estuaries embedded along the U.S. Atlantic coastline and their connections to the open ocean, such as Massachusetts Bay, Long Island Sound, Delaware Bay, and Chesapeake Bay (Fig. 1b). However, the model resolution is not fine enough to resolve smaller bays and harbors such as Tampa Bay, Galveston Bay, and New York Harbor, as well as the chains of barrier islands east of North Carolina and south of Florida to Texas. MOM6 realistically represents the broad and gently sloping continental shelf and the sharp ocean deepening across the shelf break. Previous research showed that accurate representation of coastal geometry and bathymetry is important in capturing the fine structures of storm surge (Resio and Westerink 2008; Rego and Li 2010; Mori et al. 2014).

In terms of SLR and storm surge, CM4 simulates ocean steric and dynamic effects. Like many other CMIP6 models, CM4 does not incorporate an ice sheet component, and therefore cannot simulate land ice melt and its increasing contribution to global SLR in a warming climate (Chen et al. 2017). In addition, CM4 does not simulate tides which can interact with storm surge constructively and nonlinearly and lead to the so-called tide surge (Rego and Li 2010; Muis et al. 2019). Incorporating these processes would further heighten ESL during severe storms. Other model limitations will be discussed in the discussion and conclusions section (section 6).

CMIP6 experiments with CM4 and CM4HR

As summarized in Table 1, the standard version of CM4 (1.0° atmosphere and 0.25° ocean) has been used to carry out the CMIP6 experiments including the Diagnostic, Evaluation and Characterization of Klima (DECK) and the Scenario Model Intercomparison Project (ScenarioMIP) (Eyring et al. 2016; O’Neill et al. 2016). A 250-yr model spinup was carried out prior to the DECK runs. Meanwhile, a higher-resolution version of CM4 (CM4HR) has also been configured (0.5° atmosphere and 0.25° ocean). CM4HR has been used for the High Resolution Model Intercomparison Project (HighResMIP) of CMIP6 (Haarsma et al. 2016). Daily and even hourly data critical for conducting storm and storm surge analysis have been saved. In this study, we focus on the simulations with the standard CM4 and present available results from CM4HR, thus showing the impact of atmospheric model resolution.

Table 1: CMIP6 Experiments with GFDL CM4 and CM4HR used in this study.

Daily sea level analysis

In the preindustrial control simulation (piControl) of CM4, we calculate daily mean sea level anomalies (SLA; Δhc) for a particular day and location according to

Δhc = Δηc + Δbc, (1)

where Δbc = −Δpc/ρ0g, (2)

Δηc(x,y,t)=ηc(x,y,t)−˜ηc(x,y,t1), and (3) Δpc(x,y,t)=pc(x,y,t)−˜pc(x,y,t1),t1=1,2,…,365. (4)

The subscript c denotes piControl. The terms ηc, ˜ηc, and Δηc are daily dynamic sea level (relative to a time invariant geoid), its climatology and anomaly, respectively. By definition, all of these terms have a zero global mean. Along the coast, daily fluctuations of ηc mainly reflect ocean rise and fall associated with transient weather processes and corresponding coastal waves. On interannual and longer time scales, ηc is also influenced by large-scale ocean circulation, climate modes, and external forcing. CM4 incorporates the effects on ηc of ocean temperature, salinity, and mass redistribution, as well as rainfall, evaporation, and river runoff (Griffies et al. 2014).

Because CM4 does not explicitly simulate the inverse barometer effect on sea level, we diagnose its contribution (Δbc) using sea level pressure anomalies and an equilibrium relationship [Eqs. (2) and (4)] (Ponte 2006). The terms pc,˜pc, and Δpc are daily sea level pressure and its climatology and anomaly, respectively. In this study, ˜ηc and ˜bc are removed when calculating SLA values [Eqs. (1), (3), and (4)]. But it should be noted that the absolute surge is generally higher during warm seasons than cold seasons due to the seasonal cycles (see below).

Under anthropogenic CO2 forcing, the ocean absorbs most of the excess heat due to radiative imbalance at the top of the atmosphere, and thus causing global mean thermosteric SLR (ΔGe). The term ΔGe in CM4 can be diagnosed as ΔGe(t)=−1/A∫_A ∫_{-H}^{ηe} 1/ρ0 Δρ dz dA, (5) where Δρ is the anomaly of in situ density of seawater, ρ0 is the reference seawater density, A is the global ocean surface area, H is the ocean depth, and ηe is the dynamic sea level in the CO2 experiments. The subscript e denotes CO2 experiments. In these experiments, SLA (Δhe) without global thermosteric SLR is calculated as Δhe = Δηe + Δbe, (6) where Δηe and Δbe are computed relative to ˜ηc and ˜bc in piControl: Δηe(x,y,t) = ηe(x,y,t)−˜ηc(x,y,t1), (7) Δpe(x,y,t) = pe(x,y,t)−˜pc(x,y,t1) + ϵ, t1 = 1,2,…,365. (8)

In addition to daily weather processes, regional trends of ηe and pe and the change of their seasonal cycles under the CO2 forcing also contribute to Δηe and Δbe. Note that ε is a small correction term due to the redistribution of air mass loading between the land and ocean in the CO2 experiments. SLA with global thermosteric SLR is calculated as Δhe(x, y, t) + ΔGe(t).

In summary of the analysis methods, we distinguish SLA in piControl (Δhc) and CO2 experiments (Δhe + ΔGe) to emphasize that the latter also includes global thermosteric SLR and regional dynamic SLR in addition to storm surge and other factors. Storm surge refers to the change in sea surface height relative to the predicted tide during a storm (Gregory et al. 2019). Strictly, it should not include any factors that would affect sea level in the absence of a storm. Thus, we choose to use the term “ESL” in the following to discuss high-end (extreme) daily sea levels, which incorporate all the above effects represented in CM4.

Observational and reanalysis data

In terms of data–model comparison for evaluation purposes, we use the daily tide gauge data provided by the University of Hawaii Sea Level Center (Caldwell et al. 2015) (https://uhslc.soest.hawaii.edu/). We choose long-term high-quality stations mostly along the U.S. Atlantic coast (Table 2). The data are detrended and deseasonalized. It should be noted that this comparison is not ideal since tide gauges, often located inside bays or harbors, are point measurements, while the model results represent averaged values over the coastal ocean grid cells. For the altimetry data of dynamic sea level, we use the multisatellite merged gridded dataset from the Copernicus Marine Environment Service (http://marine.copernicus.eu/). The daily anomaly data with a 0.25° spatial resolution span 1993–2017 (Pujol et al. 2016). The long-term mean dynamic sea level is based on the period of 1993–2012. For sea level pressure, we use the NCEP/NCAR reanalysis for the 1948–2018 period (Kalnay et al. 1996). The daily data have a 2.5° spatial resolution (https://www.esrl.noaa.gov).

Table 2: Daily tide gauge data used in the present study. (Note that the linear trends are not directly comparable due to different data length at different stations.)

Evaluating sea level simulations in piControl of CM4

CM4 captures the pronounced features of the long-term mean dynamic sea level observed by satellites (Figs. 3a,b). These features include the peak-to-peak range, the asymmetry associated with the gyre circulation, the sharp gradients across the Gulf Stream, Kuroshio, and the Antarctic Circumpolar Current, and the contrast between subpolar and subtropical regions and between the Pacific and Atlantic basins.

Figure 3: Dynamic sea level η in the altimetry data and piControl simulations of CM4. (a),(b) Long-term mean (m). (c),(d) Seasonal cycle as quantified by the difference between September and March (m). (e),(f) Daily variability as quantified by the standard deviation of the detrended and deseasonalized daily dynamic sea level (m). (left) Altimetry data and (right) CM4 simulations.

CM4 simulates a seasonal cycle of the dynamic sea level ˜ηc similar to the observations (Figs. 3c,d). In the Northern Hemisphere, the dynamic sea level is higher by up to 0.2 m during early autumn than during early spring, especially along the Gulf Stream and Kuroshio and nearby regions. This is mainly due to seasonal heating and cooling of the ocean, as well as seasonal changes of prevailing winds and ocean circulation. In CM4, ˜ηc along the U.S. Atlantic coast resembles that in the ocean interior, and shows increasing amplitudes from the north toward the south (Fig. 4a). In nature, annual and semiannual long tides also contribute to higher coastal sea levels during late summer and early autumn (Sweet et al. 2018). In the tropical Pacific, the belt-like feature reflects the north–south shift of the ITCZ and associated trade winds. Compared with ocean interior, ˜ηc reduces in some shelf regions such as in the Okhotsk Sea, South of Alaska, along the west coast, and on the northeast shelf of North America. The shallow ocean column on shelf limits the magnitude of seasonal thermal expansion and contraction.

Figure 4: Daily climatology and seasonal cycle of dynamic sea level ˜η and the inverse barometer effect ˜b at large coastal cities in the 150-yr piControl of CM4. (a) Dynamic sea level climatology. (b) Inverse barometer effect climatology. The long-term mean at each city is removed for better comparison. Notice the different y-axis scales in (a) and (b).

The jet-like Gulf Stream and deep western boundary current are better simulated in CM4 compared with previous model generations (Adcroft et al. 2019; Held et al. 2019). CM4 somewhat underestimates mesoscale eddy activities along the Gulf Stream, Loop Current, and Kuroshio, as well as the associated daily variability of dynamic sea level (Figs. 3e,f). This is partly due to the eddy-permitting (rather than eddy “resolving”) resolution of the ocean model. While most of the mesoscale eddies do not directly impact coastal sea levels, the warm-core rings could cause sudden TC intensification due to their anomalously high heat content (Goni et al. 2009). A notable example is Hurricane Katrina, which rapidly intensified to a category 5 hurricane after passing over a warm-core ring prior to its landfall near New Orleans (Jaimes and Shay 2009). Recent studies also suggest that in addition to direct impacts through winds and pressure, coastal storms could disrupt the Gulf Stream flow, thereby indirectly influencing sea level along the U.S. East Coast (Ezer et al. 2017; Ezer 2018, 2019).

As for surface meteorological factors, CM4 reproduces the deepening of the Icelandic low during winter and the enhanced variability of sea level pressure and surface winds along the U.S. East Coast (Figs. 5a,b). During summer, the strength and position of the North Atlantic subtropical high are realistic in the CM4 simulations (Figs. 5c,d). At higher latitudes, the summertime weather variability reduces compared with wintertime. The seasonal cycle of the inverse barometer effect (i.e., the amplitude of ˜bc) is less than 0.1 m along the U.S. Atlantic coast and its seasonal variation differs at different locations (Fig. 4b).

Figure 5: Sea level pressure (hPa) and its variability in the NCEP–NCAR reanalysis and piControl simulations of CM4. (a),(b) Mean sea level pressure (contours) and its daily variability (shading) as quantified by the standard deviation during winter (November–March). (c),(d) Mean sea level pressure and its daily variability during summer (June–October). (left) NCEP–NCAR reanalysis and (right) CM4 simulations.

Characterizing storm-related ESL in piControl

Statistics of SLA along the U.S. Atlantic coast

According to the correlation of SLA Δhc, we divide the U.S. Atlantic coast into three regions: northeast (NE), southeast (SE), and the Gulf of Mexico (GOM) (Fig. 1b). In piControl of CM4, the standard deviation σ of SLA shows a clear separation of coastal and interior ocean dynamics, roughly along the 100-m isobaths with lower σ values (Fig. 6a). Vigorous mesoscale eddies dominate in ocean interior, while wind surge and coastal waves dominate variability near the coast (Hughes et al. 2019). The bowl-shaped coastline can enhance coastal SLA variability, such as from Cape Cod to Cape Hatteras, east of Georgia, along the Florida Panhandle, and south of Louisiana and Texas. The coastal SLA variability in CM4 is consistent with the estimate from the tide gauge data, with slightly underestimated magnitudes at some sites (Fig. 6a).

Figure 6: Characteristics and statistics of SLA (Δhc) variability along the U.S. Atlantic coast in the 150-yr piControl of CM4. (a) The standard deviation (m), (b) skewness, and (c) kurtosis. Colored circles indicate the tide gauge observations (Table 2). The black line shows the 100-m isobath.

Skewness and kurtosis describe the shape of the probability distribution of SLA at different locations (see appendix A). Figure 6b shows a positive skewness of SLA along the U.S. Atlantic coast. A positive skewness indicates a longer tail at the positive end than the negative end, which occurs when ocean surge dominates in magnitude over ocean fall due to the passing of storms. TCs and ECs tend to propagate northeastward just offshore of the NE and SE coast (Fig. 2). This preferred storm track is related to the dynamics of the subtropical (Bermuda–Azores) high during summer (Elsner et al. 2000) and the midlatitude baroclinicity during winter (Lunkeit et al. 1998; Brayshaw et al. 2011). The northeasterly wind on the west and northwest side of these cyclones can cause large positive SLA values through shoreward Ekman transport. For the GOM, the northward movement (i.e., translational speed) of landfalling TCs perpendicular to the coastline causes stronger landward winds and ocean surge on the east side than the seaward winds and ocean fall on the west side. Land friction also slows down seaward winds and therefore reduces the magnitude of negative storm surge.

Kurtosis measures the “tailedness” of the SLA distribution and is a useful indicator of large storm surge and coastal vulnerability to ESL. Its value is sensitive to rare and extreme events that can lead to catastrophe. In piControl of CM4, the geographical map of kurtosis shows a striking hotspot around New Orleans with values greater than 10 (Fig. 6c). Enhanced values are also found along the west coast of Florida and the south coast of Texas. The kurtosis pattern along the coast is consistent with extreme surge height.

In addition to the moments of SLA, we also evaluate return levels of ESL by fitting the generalized Pareto distribution to the ESL values with the peak-over-threshold method (Coles et al. 2001; Arns et al. 2013) (see appendix B for details). Along the NE coast, the 1-, 10- and 100-yr return levels are generally underestimated in CM4 compared with the tide gauge data, particularly at New York (Fig. 7). The tide gauge at the Battery is inside of New York Harbor. Local large storm surge at lower Manhattan is influenced and amplified by the dynamics of the harbor and Hudson River, which are not resolved in CM4 (Blumberg et al. 2015). The ESL return levels along the SE and GOM coast are better simulated in CM4, despite the lack of the most catastrophic event like the >2-m daily surge at Galveston induced by Hurricane Carla in September 1961 (Fig. 7).

Figure 7: Comparison of SLA variability and ESL events between (left) the tide gauge data and the control simulations of (middle) CM4 and (right) CM4HR. The daily data have been detrended and deseasonalized. Six cities in the NE, SE, and GOM regions with high-quality tide gauge data are chosen: Halifax, Boston, New York (at the Battery), Charleston, Miami (at Virginia Key), and Galveston. The gray dashed lines upward show the 1-, 10- and 100-yr return levels, respectively. The shape k, scale σ, and location θ parameters (with the standard error) of the GPD fit are shown at the upper-right corner; see appendix B for details. The ESLs induced by Hurricanes Sandy and Carla are marked at New York and Galveston, respectively. Notice that the tide gauge data are point measurements at coast, while the model data are the area averaged values over the coastal ocean grid cells.

Along the GOM coast, large spikes in the simulated SLA time series are caused by landfalling TCs/hurricanes during summer and early autumn with August as the peak month (Fig. 8a), also evidenced by the long positive tail of the SLA histogram. The GOM coast is relatively quiet during winter and spring. On the NE coast, the SLA time series show a periodic wave-packet-like pattern: its variability is largely suppressed in summer but amplified in winter (Fig. 8c). Thus, most ESL events in this region occur during cold seasons associated with ECs/nor’easters (Colle et al. 2010, 2015). Nonetheless, some of the North Atlantic TCs/hurricanes can occasionally strike this northerly region during late summer and early autumn. In fact, the record high daily surge of 1.2 m at New York in the tide gauge data was induced by Sandy in October 2012 (Fig. 7), which was a large tropical–extratropical system at landfall with an unusual path perpendicular to the New Jersey shoreline (Hall and Sobel 2013). It exceeds the simulated extreme surge height at New York by CM4 and contributes to the higher 100-yr return level of ESL in the tide gauge data (Fig. 7).

Figure 8: Time series of SLA (Δhc) at New Orleans and New York in piControl of CM4. (a) SLA time series of a representative 1-yr period at New Orleans. (b) Contributions of wind surge Δηc and pressure surge Δbc to large positive surge events (Δhc > 0.2 m) at New Orleans. (c),(d) As in (a) and (b), but for New York. Note that (b) and (d) use the 150-yr piControl with seasonal cycles removed. Notice the different y-axis scales between different panels.

Wind–surge relationship

The piling up of seawater against the coast by winds is the dominant factor in storm surge. Generally, the wind effect accounts for roughly 80%–90% of the total surge height (Figs. 8b,d and 9a,b). The remainder is mainly due to the inverse barometer effect induced by the low central pressure of storms. In the following discussion, we focus on the wind surge part of the SLA.

Figure 9: Large daily surge event induced by a strong and large TC in piControl of CM4. During this event, a surge of up to 1.8 m (Δhc) occurs along the GOM coast on 24 Aug, year 138. (a) Daily dynamic sea level anomalies Δηc (m) associated with this event. (b) SLA due to the inverse barometer effect Δbc (m). (c) Daily surface wind vector and speed (m s−1). (d) Daily precipitation (cm day−1). (e) Surface current vector and speed (m s−1). (f) SST anomalies (°C) associated with the cool wake. Contours in (a), (b), (d), and (f) are daily sea level pressure anomalies (hPa) associated with the TC. The green line shows the storm track except in (c), where the line colors indicate the storm maximum daily winds (m s−1) during its propagation.

Wind stress at the ocean surface can be calculated based on the following empirical bulk aerodynamic formula: τ = Cd ρair U10 U10, (9) where τ is the wind stress vector, Cd a drag coefficient, ρair surface air density, and U10 and U10 the wind speed and vector at 10 m above the sea, respectively, computed relative to the ocean surface currents. According to Eq. (9), the wind stress magnitude is a quadratic function of 10-m wind speed. The value of Cd can increase with the increase of wind speed (Large and Pond 1981; Weisberg and Zheng 2008). At the high end of the wind spectrum associated with hurricanes and strong ECs, however, Cd reduces with the increase of wind speed (Powell et al. 2003; Oey et al. 2007).

At the coast, storm surge creates a sea surface slope and an adverse pressure gradient in the offshore (x) direction. This pressure gradient tends to balance wind stress at the ocean surface: ∂η/∂x ≈ τ/(ρ0 g H). (10) Here η is dynamic sea level, ρ0 the reference seawater density, g the gravitational acceleration, and H the ocean depth. At the coastal regions where H reduces, the sea surface slope becomes steeper and storm surge becomes higher (Pugh 1987).

The warm, semi-enclosed, and oval-shaped GOM has a basin size, geometry, and bathymetry favorable for wind setup from rotating synoptic systems (Fig. 1). When a counterclockwise rotating hurricane enters the gulf from the Caribbean Sea, water piles up at the coast due to the longshore winds and resultant shoreward Ekman transport (Hope et al. 2013). At landfall, the strongest wind is typically at its eastern and northeastern sector (Fig. 10a). So a storm track slightly west of New Orleans could realize the worst-case scenario of storm surge for the city. The intense southeasterly storm wind blows almost perpendicular into the funneling land geometry east of the Mississippi river delta, and is therefore highly efficient at raising coastal water levels (As-Salek 1998). In piControl of CM4, the wind surge height at New Orleans scales up well (r = 0.75) with the local/nearby wind speed following a quadratic relationship (Figs. 10c,d). This wind–surge relationship, as classified by the Saffir–Simpson scale (Simpson 1974; Needham and Keim 2014), highlights the nonlinear increase in coastal vulnerability as a storm intensifies.

Figure 10: Simulated wind–surge relationship in piControl of CM4. (a) Typical wind pattern for large surges at New Orleans. (b) Typical wind pattern for large surges at New York. Shading shows the correlation of large daily surge at New Orleans and New York with daily wind speed. Contours and vectors are linear regressions of sea level pressure (hPa m−1) and uυ winds (m s−1 m−1) on the large surge values. (c) Scatterplot of large daily wind surge at New Orleans and New York as a function of local wind speed. The nonlinear fit is based on a quadratic wind–surge relationship. (d) The quadratic wind–surge relationship at nine cities along the Atlantic coast. Values in the legend indicate the fit correlation.

During the landfall of a GOM hurricane, the maximum sustained wind, intense storm precipitation, and coastal surge center almost coincide (Figs. 9a,c,d). The downpour, while capable of causing inland flooding, can further increase the coastal surge height by dumping a large amount of water at the ocean surface in a short time (Wong and Toumi 2016). Over the shallow continental shelf waters, it can take a few days for the water bulge to spread and disperse through surface gravity and coastal waves. This enhancement of surge height by heavy rainfall does not work as efficiently along the NE and SE coast in part due to the narrower continental shelf. CM4 simulates the TC intensification in the GOM after passing over warm core rings north of the Loop Current (Figs. 9c,e). Due to the storm-induced ocean vertical mixing and divergent Ekman transport away from the TC center, a cool wake is evident behind the TC propagation in the CM4 simulations (Fig. 9f).

For New York, a large wind surge typically occurs when the low pressure system is located to the south and the alongshore winds induce shoreward Ekman transport (Figs. 10b,c,d). The surge at Baltimore and Miami shows weak or even no correlation with local/nearby winds. Baltimore is located in the Chesapeake Bay, where the surge is limited by the bay geometry. In nature, large storm surges in the Chesapeake Bay do exist provided that coastal storms push large amounts of waters into the bay (Ezer 2018, 2019). The narrow passage connecting the bay and open ocean is represented by one grid cell in CM4 (Fig. 1b), which may not be sufficient for simulating large inflow events. Although Miami is next to the open ocean and at the forefront of hurricane paths, the continental shelf offshore is exceptionally narrow (Fig. 1). In addition, Miami is at the southern tip of the Florida peninsula with a convex-shaped coastline. These features make storm surge less efficient at concentrating its energy. The observed tide gauge data confirm that the daily surge around Miami has not exceeded 0.4 m since the 1990s (Fig. 7). Strong winds, heavy rainfall, and big ocean waves during hurricanes are of more serious concern at Miami.

Characterizing response of ESL to CO2 forcing

Our assessment above suggests that CM4 offers a previously unavailable modeling framework to study weather–climate interactions and their combined effect on storm surge and related ESL. Next we consider a series of climate change experiments with CM4 under the CMIP6 protocol (Table 1) (Eyring et al. 2016). Among these simulations, we focus on the 1% per year increase in atmospheric CO2 concentration experiment (1pctCO2), supplemented with the companion experiments including the historical simulations, the twenty-first-century projections under two Shared Socioeconomic Pathways (SSP2–4.5 and SSP5–8.5) scenarios (O’Neill et al. 2016), and the idealized instantaneous CO2 quadrupling (abrupt 4xCO2) run. The responses of the weather–climate system and storm-related ESL are qualitatively similar between these experiments and increase in magnitude with the increase in external forcing. The results from these different experiments allow us to quantify the range of the ESL response.

Simulated changes in weather, climate, and sea level in 1pctCO2

In 1pctCO2 of CM4, both global mean surface temperature and global thermosteric SLR display upward trends during the 150-yr model simulation (Fig. 11a). Global thermosteric SLR (ΔGe) initially lags the surface temperature response, due to the gradual downward heat penetration and enormous thermal inertia of the ocean, and shows a faster acceleration after year 50. Note that ΔGe is 0.09 m at year 70 (time of CO2 doubling) and 0.34 m at year 140 (time of CO2 quadrupling); ΔGe in 1pctCO2 is corrected by removing a slow drift of the deep ocean in piControl. As a consequence of excess heat uptake mainly by the upper layers, the ocean becomes more stratified in 1pctCO2.

Figure 11: Simulated climate change and SLR in 1pctCO2 of CM4. (a) Time series of global mean surface air temperature anomalies, SST anomalies, and global thermosteric SLR. (b) Time series of dynamic SLR at large cities along the East and Gulf Coast.

In the North Atlantic and along the U.S. East Coast, ocean dynamics plays an important role in regionally modifying SLR (Levermann et al. 2005; Yin et al. 2009; Ezer 2015). In response to 1pctCO2, the Atlantic meridional overturning circulation (AMOC) weakens, which results in a 0.1-m dynamic SLR along NE at year 70 and 0.2-m dynamic SLR at year 140 in 1pctCO2 (Figs. 11b and 12a). CM4 simulates a vigorous AMOC in piControl with interannual to multidecadal variability (Figs. 13a,b). The regional enhancement of SLR along NE (on the top of global mean SLR) due to AMOC weakening is a robust feature in the previous CMIP3 (Yin et al. 2009) and CMIP5 (Yin 2012) simulations and projections, although the exact magnitude can vary across models. Compared with the previous results (Yin et al. 2009; Yin 2012; Yin and Goddard 2013), the dynamic SLR signal in CM4 extends farther southward to north of Miami (Fig. 12a). This extension may in part be due to the refined oceanic model resolution and associated representation of the continental shelf geometry and western boundary current in CM4 compared with previous model generations. Detailed mechanisms are worthy of further investigation in the future. The reduced current shear, cross-current dynamic sea level gradient, and baroclinicity tend to reduce the ocean mesoscale eddy activities (Fig. 12b).

Figure 12: Simulated ocean changes during years 131–150 in 1pctCO2 of CM4 relative to piControl. (a) Dynamic sea level anomalies (m) with a zero global mean. Contours are the long-term mean dynamic sea level (m) in piControl. (b) Response of the ocean mesoscale eddies. Shading shows the anomalies of the standard deviation of daily dynamic sea level (m). Contours are the standard deviation in piControl. To calculate the eddy-related changes in dynamic sea level, the background and large-scale SLR in 1pctCO2 is removed. (c) Pattern of SST anomalies (°C) with the global mean value removed. The green box indicates the main development region of TCs. (d) Pattern of ocean heat content anomalies (109 J m−2) with the global mean value removed.

Figure 13: Upper bound of AMOC-induced dynamic SLR under CO2 forcing. (a) Atlantic meridional overturning streamfunction (Sv) as a function of latitude and depth (m) in the long-term piControl. (b) Time series of the AMOC index defined as the maximum Atlantic overturning streamfunction values north of 30°N in piControl and the CO2 experiments. (c) Global map of dynamic sea level changes (m) during years 131–150 of abrupt 4xCO2 relative to piControl. (d) Time series of dynamic SLR in abrupt 4xCO2 at large coastal cities. The smooth curves are the exponential fit to the dynamic SLR at New York (NE), Charleston (SE), and New Orleans (GOM).

Given the importance of this dynamic SLR, it is of interest to quantify its upper bound in the stronger abrupt 4xCO2 experiment. In response to the instantaneous CO2 quadrupling, the AMOC quickly spins down and the dynamic SLR equilibrates at about 0.40 m along NE after 80 years, 0.27 m along SE, and 0.10 m along GOM (Figs. 13c,d). The e-folding time of the response is 27, 11, and 8 years, respectively. The longer response time scale at NE is likely due to the slower baroclinic processes in the higher latitudes associated with the modification of ocean density properties under CO2 forcing.

As for weather processes in a warming climate, CM4 simulates an increase in the strength (i.e., based on the maximum sustained wind and central pressure) of strong TCs/hurricanes over the North Atlantic, and a decrease in the annual count of all TCs after 100 years in 1pctCO2 (Figs. 14a,c) (Knutson et al. 2013, 2019). Despite warmer SSTs in the TC main development region (10°–25°N, 80°–20°W) (Fig. 12c), a greater warming in the tropical upper troposphere leads to a decrease in the lapse rate and an increase in static stability, which tend to inhibit atmospheric convection and TC formation in 1pctCO2 (Vecchi and Soden 2007; Knutson et al. 2008; Sobel et al. 2016). Previous studies show that the hurricane intensity may have increased during the past decades (Emanuel 2005), with stronger intensity for longer period of time (Ezer 2018). The track density map reveals that the reduction in TC frequency in CM4 mainly occurs east of the Caribbean Sea so that the number of landfalling TCs/hurricanes along the U.S. Atlantic coast remains almost unchanged (Fig. 14a). Meanwhile, extreme storm precipitation increases along the U.S. Atlantic coast in the CO2 experiments (Fig. 15), although the annual precipitation does not.

Figure 14: Response of TC and EC in 1pctCO2 of CM4. (a),(b) Changes in TC and EC track density (number per decade) during the 150-yr simulations of 1pctCO2. Shading shows the anomalies and contours show the mean in piControl. The track density is evaluated over 2° × 2° boxes. (c),(d) TC and EC count (number per year) over the North Atlantic and North America as a function of time. Notice the different scales between (a) and (b), and between (c) and (d).

Figure 15: Responses of daily (left) winds, (center) precipitation, and (right) sea level pressure anomaly along the (top) NE, (middle) SE, and (bottom) GOM coast in 1pctCO2 of CM4. The histograms use 150-yr simulations. The y axis indicates the total number of days summed over all coastal ocean grid points in the NE, SE, and GOM regions. A logarithm scale is used to better show the tail.

The total number of ECs/nor’easters offshore of the U.S. East Coast shows a more significant reduction after 100 years in 1pctCO2 of CM4 (Fig. 14b). On global scales, polar amplification of global warming can lead to a reduced meridional temperature gradient near the surface at midlatitudes, especially during wintertime (Holland and Bitz 2003; Colle et al. 2015; Shaw et al. 2016). Regionally, SST anomalies in 1pctCO2 show a “warm–cool–warm” tripolar pattern (relative to the global mean) among the regions north of the Gulf Stream, east of the subpolar gyre, and in the Nordic seas (Fig. 12c) (Rahmstorf et al. 2015). In particular, the larger ocean warming on the northeast U.S. continental shelf extends from the surface to the bottom of the shelf ocean, and is mainly caused by the weakening of the AMOC (Saba et al. 2015; Caesar et al. 2018). Recently, this region has been identified as one of the hotspots for marine heat waves (Frölicher et al. 2018; IPCC 2019) that could impact marine ecosystems (Pershing et al. 2015). This faster ocean warming, along with the faster land surface warming, reduces the temperature contrast across the Gulf Stream as well as across the land–sea boundary, thereby weakening the near-surface baroclinicity and storm track intensity of ECs/nor’easters during winter.

Compared with the SST anomalies, the maximum increase in heat content of the entire ocean column occurs at the offshore side of the shelf break (Fig. 12d). The dynamic SLR along the U.S. East Coast, the tripolar SST anomaly pattern, and the faster ocean heat accumulation along the shelf break are consistent manifestations and consequences of the AMOC weakening in CM4.

Response of ESL to CO2 forcing

Figure 16 compares the SLA distribution in the historical, SSP projection, 1pctCO2, and abrupt 4xCO2 runs with (Δhe + ΔGe) and without (Δhe) global thermosteric SLR. It is evident that the increase in CO2 forcing progressively shifts the probability density function (PDF) curve to the right and toward higher values. In the historical run, the shift is relatively small due to the anthropogenic aerosol forcing largely counteracting greenhouse gas forcing until about 1990, which leaves insufficient time for sea level response (Held et al. 2019). The shift is more significant in the historical run without anthropogenic aerosol and land use forcing, in the SSP projections, and in 1pctCO2, and is strongest in the abrupt 4xCO2 run.

Figure 16: Progressive responses of ESL to different strength of external forcing. The histograms indicate the distributions of SLA values in piControl, 1pctCO2, abrupt 4xCO2, the historical runs (with and without anthropogenic aerosol and land use forcing), and SSP projection runs of CM4. (a) SLA in piControl (Δhc) and climate change experiments without global thermosteric SLR (Δhe). (b) SLA with global thermosteric SLR (Δhe + ΔGe). The y axis indicates the total number of surge days summed over all coastal ocean grid points in the NE, SE, and GOM regions. All values have been normalized to a 150-yr period for comparison. (left) The piControl with idealized 1pctCO2 and abrupt 4xCO2 simulations and (right) piControl with the historical and SSP projection runs, for (top) NE, (middle) SE, and (bottom) GOM.

Without global thermosteric SLR, the nearly uniform shift of the PDF curve in 1pctCO2 and the elevated daily surge height along the NE and SE coast is mainly attributable to the AMOC-induced dynamic SLR (Fig. 16a). The Kolmogorov–Smirnov statistical test indicates that the shift of the PDF curve is statistically significant. Along the GOM coast, the overall shift of the PDF curve to the right is relatively small, except having a disproportionately longer tail (Fig. 16a). This heightening of ESL is consistent with the increase in TC intensity under the CO2 forcing (Fig. 15). Adding global thermosteric SLR substantially widens the SLA distribution, and reduces its skewness and kurtosis (Fig. 16b).

In piControl of CM4, the return levels for 1-, 10- and 100-yr ESL events differ dramatically among three major coastal cities (Fig. 17) (see appendix B). The tightly packed return levels at Miami are lowest, in sharp contrast with the highest and widely separated return levels at New Orleans, especially for the 1-in-100-year events (0.26 m at Miami vs 1.83 m at New Orleans) (Tebaldi et al. 2012). SLA at Miami shows small variability and a lack of tail at both ends of its histogram (Fig. 17b). The opposite occurs at New Orleans with large surge spikes and a long histogram tail (Fig. 17c), while the surge at New York is in the middle (Fig. 17a).

Figure 17: Time of emergence of the anthropogenic signal in ESL in 1pctCO2 of CM4, for (a) New York, (b) Miami, and (c) New Orleans. Blue and red colors denote SLA in piControl (Δhc) and 1pctCO2 (Δhe + ΔGe), respectively. Horizontal dashed lines denote the 1-, 10-, and 100-yr return levels of daily sea level in the 150-yr piControl based on the GPD method. Triangles and diamonds indicate TOE in ESL height and frequency, respectively. Rectangles denote permanent exceedance by the rising mean sea level. Shown are the (left) time series and (right) histogram of SLA.

In climate change studies, the time of emergence (TOE) of the anthropogenic signal is an important quantity for detection and attribution purposes (Diffenbaugh and Scherer 2011; Hawkins and Sutton 2012). Under CO2 forcing, the anthropogenic signal can emerge in terms of ESL height or frequency or both. With 1pctCO2 of CM4, we quantify and compare TOE in terms of ESL height and frequency with and without global thermosteric SLR (see appendix C for the TOE calculation method). With global thermosteric SLR (Δhe + ΔGe), TOE in ESL height of the 1-yr events occurs at year 23, 22, and 70 for New York, Miami, and New Orleans, respectively (Fig. 17). It is longer and later for the 10-yr events, and occurs at year 69 and 50 for New York and Miami, respectively. At New Orleans, the 10-yr signal emerges in ESL frequency (at year 86) rather than in ESL height. For the more extreme 100-yr event, TOE in ESL frequency can be identified at year 64 and 55 for New York and Miami, respectively. However, the 100-yr signal cannot be detected at New Orleans. This is mainly due to the large natural variability and a slower SLR at New Orleans due to ocean steric and dynamic effects in CM4.

The early TOE in ESL at New York is facilitated by the AMOC-induced dynamic SLR in this region (Figs. 11b and 12a). The anthropogenic signal first shows up in wintertime ESL events and coastal flooding associated with nor’easters. The early TOE at Miami is primarily due to the weak background SLA variability especially the low surge height in piControl. More significantly in 1pctCO2, the 1-, 10- and 100-yr return levels of ESL are permanently exceeded at Miami by the rising mean sea level at year 71, 92, and 102, respectively (Fig. 17b). By contrast, New Orleans shows no permanent exceedance. We find that along the U.S. East Coast, the anthropogenic emergence occurs even before global thermosteric SLR is taken into account.

Impact of the atmospheric model resolution

The refined atmospheric model resolution (0.5°) in CM4HR leads to more intense TCs/hurricanes with stronger winds and smaller size in the control simulation (Fig. 18). These magnify SLA extremes at both ends, but with a greater influence on the positive side along the SE and GOM coast (Fig. 19). For example, the highest daily surge at the GOM coast increases from 1.8 m in CM4 to 2.3 m in CM4HR. The simulations of ECs in CM4 and CM4HR are similar due to their large size relative to TCs. Compared with CM4, the return levels of the 1-, 10- and 100-yr ESL are higher in CM4HR, closer to those for the tide gauge data (Fig. 7). We find that in abrupt 4xCO2, the responses of storm-related ESL and the impact factors are qualitatively similar between CM4 and CM4HR, including storm characteristics, surge statistics, and oceanic and atmospheric circulation, as well as global and regional sea level (Fig. 19).

Figure 18: Large daily surge event induced by a strong TC in the control run of CM4HR. During this event, a surge of up to 2.3 m (Δhc) occurs along the GOM coast on 1 Sep, year 80. (a) Daily dynamic sea level anomalies ηc (m) associated with this event. (b) SLA due to the inverse barometer effect Δbc (m). (c) Daily surface wind vector and speed (m s−1). (d) Daily precipitation (cm day−1). (e) Surface current vector and speed (m s−1). (f) SST anomalies (°C) associated with the cool wake. Contours in (a), (b), (d), and (f) are daily sea level pressure anomalies (hPa) associated with the TC. The green line shows the storm track except in (c), where the line colors indicate the storm maximum daily winds (m s−1) during its propagation. Note that Fig. 9 shows similar plots but for CM4.

Figure 19: Histograms of SLA in the 150-yr control and abrupt 4xCO2 runs of CM4 and CM4HR. (left) SLA in the control run (Δhc) and in the CO2 run without global thermosteric SLR (Δhe) and (right) SLA in the control run (Δhc) and in the CO2 run with global thermosteric SLR (Δhe + ΔGe), for (top) NE, (middle) SE, and (bottom) GOM.

Conclusions and model limitations

In the present study, we use a seamless and self-consistent global modeling framework (GFDL CM4) to study weather–climate interactions and their combined effect on extreme sea level along the U.S. Atlantic coast. Thanks to recent progress in model development and improvement, some outstanding questions of significant societal concerns can be answered now for the first time. We compare the characteristics of storm-related ESL among the NE, SE, and GOM regions and their responses to CO2 forcing. We find that under internal weather processes, the low-lying Gulf Coast is most vulnerable to hurricanes and related storm surge. New Orleans is a striking hotspot with the highest surge efficiency in response to storm winds. In response to a 1% per year atmospheric CO2 increase, the elevated surge height along the U.S. East Coast is mainly caused by the background SLR, while that along the Gulf Coast is sensitive to the modification of hurricane characteristics by the external forcing.

Our results confirm previous findings (e.g., Yin et al. 2009) that among the densely populated coastal regions worldwide, the U.S. East Coast is special and more vulnerable in terms of dynamic SLR (Fig. 13c). The AMOC-induced regional SLR facilitates the early emergence of the anthropogenic signal in daily surge, especially during wintertime flooding associated with nor’easters. The weakening of AMOC in the CO2 experiments is mainly caused by thermohaline processes (Gregory et al. 2005; Stouffer et al. 2006; Hu et al. 2009; IPCC 2019). On shorter time scales, recent research showed that it is the atmospheric wind and pressure that influenced annual mean sea level along the U.S. northeast coast (Piecuch et al. 2019). While different possible factors need to be explored, our results here from the new CMIP6 simulations stress that the active and important role of AMOC in weather, climate, regional dynamic SLR, and storm-related ESL should not be underestimated, particularly for the twenty-first century.

Nonetheless, given the complexity of SLR and storm surge along the U.S. Atlantic coast, there are important caveats about model limitations. In nature, the highest water level typically occurs during tide surge. In addition to tidal ranges, the surge–tide nonlinear interactions depend on multiple factors such as the timing of landfall, the distance to the storm, and the slope of the continental shelf (Rego and Li 2010). CM4 does not simulate tides as well as wave setup or run-up, therefore underestimating the highest water level during storm surge. CM4 does not implement a wetting and drying scheme to represent the intrusion of seawaters and coastal inundation during storm surge (Hubbert and McInnes 1999). The ESL analysis based on the daily mean data can underestimate the peak hourly surge. Uncertainties also come from the lack of the strongest (e.g., category 4 and 5) hurricanes in CM4 and CM4HR, as well as the underestimated return levels of ESL (Fig. 7). An increase in the return level in piControl could delay TOE in 1pctCO2.

Without an ice sheet model, CM4 cannot simulate the impacts of Greenland melt on sea level, AMOC and geoid changes along the U.S. East Coast (Kopp et al. 2010; Bakker et al. 2016). Finally, CM4 does not simulate climate-unrelated factors. Most of the U.S. Atlantic coast is influenced by land subsidence, particularly at New Orleans and along the Texas coast (Nienhuis et al. 2017), which can increase relative SLR and exacerbate the impact of storm surge and coastal flooding (Table 2). By contrast, land uplift in some of the New England coastal regions can offset and mitigate the dynamic SLR (Karegar et al. 2016).

Ideally, projections of SLR and storm-related ESL along the U.S. Atlantic coast should take all these factors into account. We trust that future model development will continue to address these and other limitations and further improve the model’s ability and capacity, thereby providing more accurate information for effective planning and preparedness along the U.S. Atlantic coast.

Appendix

Appendix A: Calculation of Skewness and Kurtosis of SLA

Skewness and kurtosis describe the shape of the probability distribution of SLA (Hughes et al. 2010). In piControl, the skewness of SLA values is the third standardized moment

skew = E [ ( Δhc − μ / σ )^3 ], (A1)

with μ and σ being the long-term mean and standard deviation of Δhc, and E the expectation operator. The kurtosis is the fourth standardized moment and is computed by

kurt = E [ ( Δhc − μ / σ)^4 ]. (A2)

The kurtosis of a normal (Gaussian) distribution is 3. Sometimes it is useful to compare the kurtosis of a distribution to this value (the so-called excess kurtosis). In this study, we use the original calculation of kurtosis based on Eq. (A2). Large kurtosis values indicate more extreme outliners in the distribution.

Appendix B: Return Level and Period of ESL in piControl

We use the peak-over-threshold method (Coles et al. 2001; Arns et al. 2013) to calculate return levels of storm-related ESL corresponding to particular return periods in the 150-yr piControl. We set the 99th percentile of SLA (Δhc) as the threshold to extract the subset of extreme values. We fit the empirical distribution of the subset with the generalized Pareto distribution (GPD).

y = p(x|k,σ,θ): = (1/σ) (1 + k (x − θ)/ σ)^{−1−1/k} , k ≠ 0 = (1/σ) e^{−(x−θ)/σ}, k = 0, (B1)

for x ≥ θ when k > 0, and θ ≤ x ≤ θ − σ/k when k < 0. Here y is the PDF of GPD; k, σ, and θ are the shape, scale, and location parameters, respectively. For k = 0, GPD becomes the exponential distribution. Given a particular return period T (1, 10, or 100 years) of ESL events, the corresponding return level zT is

zT = θ + σ/k { [ (1 − F(θ)) / p ]^k − 1 }, k ≠ 0: = θ + σlog [ (1 − F(θ)) / p ], k = 0 (B2)

where p and F(θ) are the probability of the event and cumulative density function, respectively.

Appendix C: Time of Emergence of ESL in 1pctCO2

For the 1-in-1-year ESL events, we use a 20-yr moving window and the 99.73 percentile of SLA (empirical 1-yr event return level) as the threshold to extract the event samples (Si, i = year 20, 21, …, 150). For example, S20 and S21 denote the subset exceeding the 99.73 percentile of SLA during years 1–20 and 2–21, respectively. We then identify the median value of Si in piControl (m^i_c) and 1pctCO2 (m^i_e). The subscripts c and e denote piControl and CO2 experiments, respectively. The mean Mc and standard deviation σc of m^i_c are

Mc = 1/131 ∑^150_{i=20} m^i_c,

σc = [∑^150_{i=20} (m^i_c - Mc)^2 / 131 ]^{1/2}. (C1)

In 1pctCO2, TOE in ESL height of the 1-yr events is defined as the year i beyond which m^i_e permanently exceeds (Mc + 2σc). If the anthropogenic signal emerges within the first 20 years, TOE is marked as year 20. For the 10-yr events, we follow the same detection procedure except using a 50-yr moving window instead to increase the sample size and the 99.97 percentile as the threshold. We do not evaluate TOE in ESL height for the 100-yr event due to its extreme rareness in the 150-yr simulations of CM4.

To quantify TOE in ESL frequency, we count the occurrence time above the 1- and 10-yr return level in piControl using a 20- and 50-yr moving window, respectively. TOE in ESL frequency is defined as the year beyond which the occurrence time in 1pctCO2 permanently exceeds the mean plus two standard deviation of the occurrence time in piControl (similar to the method above for ESL height). For the 100-yr event, TOE is the year beyond which its occurrence time in 1pctCO2 within a 50-yr moving window permanently exceeds 1 (i.e., ≥2).