Projected Changes to Cool-Season Storm Tides in the 21st Century along the Northeastern United States Coast

,

While often less severe than hurricane-driven surge, ETCs are responsible for most moderate surge events in the NEC region and affect a very wide region of the coastline (Booth et al., 2016). It is thus important to assess the impacts of ETC-driven surge and storm tides, especially when considering events of moderate frequency (1-years to 3-years time scales). Systematic changes to ETCs under a changing climate could affect the frequency and severity of cool-season storm tides leading to more (or less) frequent coastal flooding. For instance, it has been shown that large ETC-driven surges along the NEC are typically generated by slow-moving deep cyclones to the south of a strong anticyclone (Catalano & Broccoli, 2018), so changes in these types of events would be play a critical role in altering the frequency of ETC-driven coastal flooding events. Furthermore, any changes to storm tides must be viewed with respect to rising sea levels that would further enhance the risk to coastal flooding (Booth et al., 2016), so it is also important to put the magnitudes of each into context.
Several studies have explored the effects of global warming on cool-season ETC climatology for the NEC. Most of these project a reduction in the density of ETCs over the continental United States and western North Atlantic Ocean (Chang, 2013;Colle et al., 2013;Long et al., 2009;Seiler et al., 2018;Teng et al., 2008). ETC intensities over the western North Atlantic Ocean are also predicted to weaken; however, cyclones may become more intense and deepen more rapidly just inland of the NEC (Colle et al., 2013). Many of these studies used output from global climate models (GCMs), particularly those from Phase 5 of the Coupled Model Intercomparison Project (CMIP5; Taylor et al., 2012) which have horizontal model resolutions of ∼100-300 km. Due to the regional differences and dependence on model resolution, Colle et al. (2013) suggests that dynamically downscaled regional versions of the GCMs are needed to investigate the changes to ETC track density and intensification in more detail. Other studies show that 20-km horizontal resolution dynamically downscaled regional climate models (RCMs) generate stronger cyclones than the parent GCMs based on the surface wind speed (Booth et al., 2018;Zhang & Colle, 2018). The downscaled simulations also indicated that latent heating related diabatic processes, which are otherwise too weak in coarse-resolution GCMs, could enhance development of intense ETCs over the NEC (Zhang & Colle, 2018). However, studies by Long et al. (2009) and Seiler et al. (2018) suggest that projected changes to ETC density are not particularly sensitive to model resolution.
Previous studies have examined climate change impacts on cool-season surge along the NEC using statistical (Roberts et al., 2017) or hydrodynamic  surge models forced by surface winds and pressure from CMIP5 GCM ensembles. Roberts et al. (2017) found no significant change to surge return intervals in a future period (2054-2079) compared to a historical period  at The Battery in New York City. This was attributed to the fact that projected ETC changes did not occur in regions that favor the generation of surge at The Battery. Similarly, Lin et al. (2019) found relatively small projected changes (<7%) to extreme storm surge heights for the same future period along most of the NEC, while noting however that one of the GCMs showed a more substantial increase of up to 36% for the 50-years surge height. These previous analyses contain uncertainties due to the usage of atmospheric forcings from coarsely resolved GCMs and surge prediction by straightforward multilinear regression (Roberts et al., , 2017 or a relatively coarse-resolution hydrodynamic model (∼1-km coastal resolution; Lin et al., 2019). Furthermore, astronomical tides were omitted even though it is the combination of surge and tide (storm tide) that needs to be considered to assess local flooding potential (Horsburgh & Wilson, 2007). In this study, we address these limitations by integrating results from three 12-km dynamically downscaled RCMs with a high-fidelity (∼50-m coastal resolution) hydrodynamic storm tide model, which we run multiple times to account for the arbitrary surge-tide phasing. Using this high-resolution integrated modeling system, we aim to: (1) quantify projected changes and associated uncertainties of cool-season storm tides in the 21st century along the NEC as compared to estimates of sea level rise (SLR) and (2) relate projected storm tide changes to the ETC climatology, such as changes to track patterns and intensity.
The rest of this paper is organized as follows: first, we introduce the modeling and analysis approach in Section 2; in Section 3, we present results showing model accuracy during the historical period, followed by our projected changes to storm tides in future decades and the associated ETC patterns driving these storm tide changes; in Section 4, we discuss our major findings and their implications, as well as the uncertainties and limitations of the study.

Dynamically Downscaled Regional Climate Model Experiments
Three sets of CMIP5 GCMs (CCSM4, GFDL-ESM2G, and HadGEM2-ES) have been dynamically downscaled to 12 km horizontal resolution using the Weather Research and Forecasting (WRF) v3.3.1 model (Wang & Kotamarthi, 2015;Zobel et al., 2018). These three GCMs were chosen based on evidence that they approximately represent the spread of climate sensitivity for the 30 GCMs in the CMIP5 experiment (GFDL-ESM2G-lower sensitivity, CCSM4-moderate sensitivity, HadGEM2-ES-high sensitivity) (Sherwood et al., 2014;Zobel et al., 2018). Herein, the downscaled RCMs are referred to as WRF-CCSM4, WRF-GFDL, and WRF-HadGEM. The WRF computational domains cover all of the continental USA (see Figure S1) extending out to the western North Atlantic Ocean to encompass all of the U.S. East Coast and Gulf Coast, and parts of the Caribbean (Figure 1). The mean of certain atmospheric variables (zonal and meridional wind, geopotential height, temperature, and relative humidity) in the CCSM4 and GFDL-ESM2G GCMs were bias-corrected before using these as boundary conditions for the RCM, while the HadGEM2-ES GCM was not bias-corrected. Spectral nudging at the lateral boundary was also applied to the WRF-CCSM4 RCM. The impacts of the GCM bias-correction and spectral nudging were explored in Zobel et al. (2018). Sea level pressures (SLP) and 10-m wind velocities (U10; both zonal and meridional directions) are output from the model simulations at 3-hourly intervals and used to force the hydrodynamic storm tide model.
Each RCM provides meteorological data for three decadal periods; 1995-2004 ("historical" decade), 2045-2054 ("midcentury" decade), and 2085-2094 ("late-century" decade). This corresponds to nine continuous cool-seasons for each decade, e.g., November 1995to March 1996to November 2003to March 2004 for the historical decade. The RCM data set is limited to decadal periods due to the large computational requirements of the high-resolution WRF model that encompasses the entire continental USA. The future decades were simulated under the Representative Concentration Pathway (RCP) 8.5, a pathway that assumes high levels of greenhouse gas emissions by 2100 with an effective radiative forcing increase of 8.5 W/m 2 due to a large global population and little technological improvement (Riahi et al., 2011). Recent estimates predict that RCP8.5 will accurately represent current emissions out until midcentury and represents at least plausible levels out until late-century (Schwalm et al., 2020).

Hydrodynamic Storm Tide Model
The ADvanced CIRCulation (ADCIRC) hydrodynamic model (Luettich & Westerink, 2004) is used to simulate storm tides along the NEC. We use Version 55 of the model that newly incorporates self-attraction and loading (SAL) tides, internal tide induced wave drag, and modifications to the governing equations to correctly account for Earth's curvature (Pringle et al., 2021). The ADCIRC computational domain covers the western North Atlantic Ocean west of the 60° meridian (Figure 1a), a well-studied region for the ADCIRC model (e.g., Bunya et al., 2010;Hope et al., 2013;Marsooli & Lin, 2018;Roberts et al., 2019b;Westerink et al., 2008). Version 3.0.0 (Pringle & Roberts, 2020) of OceanMesh2D (Roberts et al., 2019a) is used to automatically generate an unstructured mesh for the study domain using carefully designed combinations of shoreline geometry and seabed topography-based element sizing functions (cf. Roberts et al., 2019b). A nominal minimum element size concentrated at the coast is set to 50 m in the NEC region and 1 km elsewhere ( Figure 1). The nominal maximum element size in the deep ocean is set to 10 km. Mesh bathymetry is interpolated from the high-resolution (∼1-3 m) USGS Coastal National Elevation Database (CoNED) in the NEC region and ∼500 m SRTM15+ (Tozer et al., 2019) Version 2 data elsewhere. The storm tide model is forced with SLP and U10 (both zonal and meridional directions) from the downscaled WRF climate model data, in addition to astronomical tidal potential and SAL for the eight dominant tidal constituents (M 2 , S 2 , N 2 , K 2 , K 1 , O 1 , P 1 , Q 1 ). Astronomical tides are also prescribed at the open boundary using the TPXO9-Atlas (Egbert & Erofeeva, 2019). To account for the random timing between tides and the surge-producing weather event we simulate each season five times with different tidal phases (−10, −5, +0, +5, +10 h offsets from the actual date-time). A computational time step of 12 s was used for all simulations, and water elevations were output at 1-h intervals for the analysis.

Peak Storm Tide Elevations
For each decade and each realization of the five tidal phases, we extracted Peak Storm Tide elevations (PST) separated by a minimum of 3 days from the data to identify unique ETC-driven events . In previous studies, this data has been processed into extreme value estimates of low frequency PST events, e.g., 50-years and 100-years return periods, obtained by fitting the tail of extracted peaks to the Generalized Pareto Distribution using the Peak Over Threshold method Marsooli et al., 2019). However, we deemed the decadal-long simulations in this study to be too short to conduct a robust extreme value analysis. We instead choose to measure changes in PST for return periods contained within the time period of the simulations; the 3-season and 1-season return periods. We define the 3-season PST empirically as the third highest PST within a decade (i.e., the third largest in nine cool seasons); while the 1-season PST is defined as the ninth highest PST within a decade.
Simulated PST values are reduced to a single value for each county along the NEC coast ( Figure 1) so that the results are more easily presented and understood (these results are described in Sections 3.2 and 3.3). The value for each county is taken as the maximum PST at the mesh vertices along the coastline of that county (cf. Marsooli et al., 2019). For comparison, the astronomical MHHW (mean higher high water) value for each county is also approximated from harmonic constituent amplitudes (≈1.1M 2 + K 1 + O 1 -half of the sum of the mean range and diurnal range, Parker, 2007). Differences in the county-wide PST values between future and historical decades are presented individually for each RCM forcing. The tidal phase related uncertainty in the difference is found by taking the minimum, mean, and maximum differences of all possible combinations of tidal phase in the future and historical decade (25 total). Furthermore, we compare the relative magnitude of storm tide changes to SLR under the RCP8.5 scenario, which is computed for each county from the ocean model outputs of the three parent GCMs. SLR is approximated as the difference between the future cool-season decadal average and the historical cool-season decadal average of the total sea surface height in the closest GCM ocean point to the county midpoint. We define the total sea surface height as the sea surface height (CMIP5 variable zos) plus the global average steric sea level change (CMIP5 variable zossga) (Becker et al., 2016)

Cyclone Tracking and Mapping to Peak Storm Tides
To attribute changes in storm tides to patterns of ETC tracks and intensities, we extracted storms from the meteorological data by tracking the local minimums of SLP using Version 2 of CycloneTrack (Flaounas et al., 2014), a cyclone tracking algorithm. To filter out small scales in SLP, a 2-D Gaussian smoothing kernel with a standard deviation of 10 is used.
We select ETC tracks that produce the nine highest peak storm tide elevations for each RCM within one of the following six multicounty subregions: New England (NE), Long Island Sound (LIS), New York/New Jersey Bight (NY/NJ), Delaware Bay (DB), Chesapeake Bay (CB), and North Carolina (NC) (Figure 1). Tracks are selected by finding those that exist within the NEC domain just before and after the time of the peak storm tide. Usually there is just one track that meets these criteria, but sometimes there are no tracks in which we skip to the next highest peak storm tide elevation, or very occasionally there are two tracks in which we record both. The results of this analysis are presented in Section 3.4.

Historical Data
0.25° ERA5 reanalysis data (European Centre for Medium-Range Weather Forecasts, 2019) are used to evaluate the downscaled RCMs during the historical decade. Note that the ERA5 reanalysis data were sampled at 3-hourly intervals to match the RCM model outputs. Comparisons are made based on the statistics of wind velocities and ETC frequency and spatial distribution (see Section 3.1.1). Furthermore, ERA5 was used to force the hydrodynamic in the same way that it was forced by the downscaled RCMs for comparison (Section 3.1.2).
Observed water levels at National Oceanic and Atmospheric Administration (NOAA) tide gauges are used to evaluate the accuracy of the simulated storm tides by the hydrodynamic model via errors in the daily maximum water levels (DMWLs) and PSTs under ERA5 and downscaled RCM forcing (see Section 3.1.2). Note that the water elevations in the hydrodynamic model are "storm tide elevations" which we mean here to be the barotropic response to astronomical and atmospheric forcing and are not referenced to a particular vertical datum as sea level rise and other low frequency components to sea level are missing. Therefore, to compare the observed tide gauge data to the model, the mean of the observed water elevations for each season is subtracted from the observed water levels to approximately obtain the observed "storm tide elevation."

Dynamically Downscaled Regional Climate Model
Figures showing the historical accuracy of the WRF-based RCM simulations compared to ocean buoy observations and ERA5 reanalysis data are presented in the supplementary material. Low-level winds in the RCMs during the historical decade are shown to be mostly accurate (RMSE < 0.6 m/s) at offshore buoy locations and in the northern part of the NEC region, while errors are largest (RMSE up to 2.4 m/s) near the North Carolina and Virginia coastline ( Figure S2). Higher errors at the coast could be due to land-masking and geopotential height discrepancies in the WRF model. WRF-HadGEM, which does not have GCM bias-correction, has greater wind speed errors than the other two RCMs at offshore locations, but WRF-HadGEM errors at buoys closer to the coast are in fact smaller at some locations. Figures S3-S5 show median and 95% quantile wind speeds and directions for the RCMs. At many of the buoy locations all three RCMs tend to overestimate winds speeds coming from the south, while there is generally better agreement for easterly, westerly, and northerly winds.
Compared to ERA5, the density distribution of simulated ETCs in the NEC region is shown to be accurate to within 1-2 cyclones/season for all the RCMs ( Figure S6), showing improvements over the parent GCMs which tend to underestimate cyclone densities further (Roberts, 2015). WRF-CCSM4 underestimates the density of ETCs offshore while WRF-HadGEM overestimates. WRF-GFDL overestimates cyclone density closer to the coast especially near the New Jersey and New York region. Compared to ERA5, WRF-CCSM4 shows good agreement in the shape of the distribution of ETC maximum lifecycle intensities [minimum SLP (P min ) and maximum U10 (U max )], but it produces too few of the most frequently observed cyclones (P min ∼ 990-1,010 hPa, U max ∼ 14-24 m/s; Figure S7). The WRF-GFDL and WRF-HadGEM produce a greater overall number of ETCs than WRF-CCSM4, which matches more closely with ERA5 ( Figure S7). However, there are more ETCs of greater intensity (P min < 990, U max > 24 m/s) than in ERA5. Similarly, Figures S3-S5 show that the 95% quantile wind speeds against ocean buoys are more overestimated in WRF-GFDL and WRF-HadGEM RCMs than in WRF-CCSM4.

Hydrodynamic Storm Tide Model
First, root-mean-square errors (RMSE) from quantile-quantile plots of DMWLs at NOAA tide gauges are presented to demonstrate the historical accuracy of the hydrodynamic for everyday conditions. Errors are within 0.2 m throughout most of the NEC ( Figure 2); the tidal range at one tide gauge in a complex wetland environment (Wachapreague, VA) is overestimated by the model leading to the highest errors there (RMSE ∼ 0.3 m). The errors of the RCM-forced runs are distributed similarly to those under the ERA5 reanalysis forcing, partly because both surge and tidal components contribute to the DMWL. The greatest difference in DMWL errors between the atmospheric forcings is found in Delaware Bay and Chesapeake Bay. Wind speed errors ( Figure S2) just offshore of Delaware Bay and Chesapeake Bay were shown to be largest for the WRF-HadGEM model which could be contributing to the also generally larger associated DMWL errors in these estuaries. Wind speed errors inside Chesapeake Bay is greatest for WRF-CCSM4 and associated DMWL errors throughout Chesapeake Bay are indeed larger than those under ERA5 and WRF-GFDL forcing.
The simulated 3-season PST (mean of the five tidal phase realizations) during the historical decade at each mesh vertex in the NEC region under ERA5 reanalysis and RCM forcing is shown in Figure 3. Furthermore, the 3-season and 1-season PST are compared to selected NOAA tide gauges in the NEC region in Figure 4 using error bars to indicate the tidal phase-based un certainty, whereas the observational data (demeaned for each season) is a single realization of the tidal phase. This uncertainty is greater for the lower frequency 3-season PST (Figure 4a) than the 1-season one (Figure 4b), and is greatest along the stretch of coastline from Woods Hole to Atlantic City tide gauges under WRF-GFDL forcing. Although the spatial distribution of PSTs is similar under all meteorological forcing since it is heavily related to the tidal range, the RCM forcing consistently generates greater storm tide elevations than the reanalysis forcing. The RCM-forced PSTs are indeed overestimated at some tide gauges, in particular, those in the Long Island Sound and Delaware Bay (Woods Hole to Montauk, Philadelphia, and Reedy Point). However, it is equally true that PSTs under reanalysis forcing are underestimated at many of the tide gauges, which could mean that the most extreme winds are smoothed-out in the reanalysis data compared to the RCMs.

Peak Storm Tide Elevation Changes in Midcentury Decade
Projected PST changes range within ±0.4 m for the 3-season PST (Figure 5b) and ±0.3 m for the 1-season PST (Figure 6b) throughout the NEC region for the three RCM forcings by midcentury. The magnitude of the 3-season PST changes by midcentury across the NEC is of similar order to SLR under the CCSM4 and GFDL SLR scenarios (∼0.2 m). However, the magnitude of SLR according to HadGEM (0.6-0.8 m) is significantly larger than 3-season PST changes. The magnitude of 1-season PST changes by midcentury across the NEC is slightly smaller in magnitude to CCSM4/GFDL SLR (∼0.2 m), and much smaller than HadGEM SLR (0.6-0.8 m).
In the northern counties (<#60), projected PST changes are highly dependent on the RCM forcing. WRF-CCSM4 forcing results in mostly small average (tidal phase-based) increases to PSTs in the New England, Long Island Sound, and New York/New Jersey Bight subregions and up to as large as 0.15 m along the Delaware River (#47-55) in the Delaware Bay subregion. Similarly, WRF-HadGEM forcing projects on average mostly  For southern counties (>#60), WRF-CCSM4 and WRF-HadGEM forcings largely project average decreases. While little change to PSTs is seen in the central west regions of Chesapeake Bay (counties #60-90), fairly large average decreases are projected under these two RCM forcings elsewhere, up to 0.2-0.3 m. Comparatively, the WRF-GDFL run projects mostly small average increases to PSTs of ∼0.2 m for the 3-season and ∼0.1 m for the 1-season in these southern counties. Localized larger average increases of ∼0.15-0.2 m to the 1-season PST are noticeable at counties located up western tributaries of Chesapeake Bay and Albemarle Sound (#86-88, #94-95, #103, #121-122). The tidal phase-based variability is not very significant for these southern counties, especially for those in the north part of the North Carolina subregion. This area corresponds to Albemarle and Pamlico Sounds that have a tidal range of only a few centimeters (refer MHHW in Figure 5a), and hence susceptible to changes in PST due to locally generated storm surge only.

Peak Storm Tide Elevation Changes in Late-Century Decade
Projected PST changes range within ±0.5 m for the 3-season PST (Figure 5c) and ±0.3 m for the 1-season PST (Figure 6c) throughout the NEC region for the three RCM forcings by late-century. In this decade, the WRF-GFDL and WRF-HadGEM forcing project changes that have a high-degree of spatial variability, while WRF-CCSM4 shows smaller and smoother changes. The magnitude of the 3-season PST changes in late-century across the NEC is similar to but mostly smaller in magnitude to SLR under the CCSM4 and GFDL SLR scenarios (∼0.4 m). The magnitude of SLR according to HadGEM (1.1-1.3 m) is far greater than the 3-season PST changes; it is at least two times the size of the greatest PST decrease for any run and county. The magnitude of the

Cool-Season Storm Climatology Patterns Driving Storm Tide Changes
In the northern subregions (New England, Long Island Sound, New York/New Jersey Bight) during the historical period, many of the RCM-simulated highest nine peak storm tide generating ETC tracks tend to follow nearby and parallel to the coastline (Figure 7). In the future decades, ETC tracks tend to be more sparsely distributed and are less likely to follow that coastline parallel track, either veering further offshore or tracking more inland in the south-to-north direction. Specifically, in WRF-GFDL, there are fewer storms that make close or direct impact to the subregion of concern which could explain why WRF-GFDL forced PSTs are projected to mostly decrease in the northern subregions in future decades. However, WRF-GFDL forced PSTs are projected to increase in the Hudson River in the late-century decade. It is noticeable that all WRF-GFDL decades produce storms that pass Figure 7. Distribution of regional climate model (RCM)-simulated cool-season ETC tracks that produce the nine highest peak storm tide elevations under each RCM forcing within a northern subregion as indicated by the black dashed-line polygon (top: New England, middle: Long Island Sound, bottom: New York/New Jersey Bright). The tracks are color coded according to the RCM and line thicknesses are proportional to ETC intensity (P min , lower is thicker). Tracks are shown for the three decades investigated in this study (left-to-right): historical (1995-2004), midcentury (2045-2054), and late-century (2085-2094). ETCs, extratropical cyclones. through the northern New York/New Jersey Bight subregion where Hudson River is located, but in the late-century, the storm tracks are clustered further to the north than the other two decades. Due to the anticlockwise cyclone rotation, southerly winds which would be favorable to surge generation in the river are to the southeast of the cyclone center. Similarly, WRF-HadGEM shows a number of relatively intense storms tracking just inland of the Long Island Sound and New York/New Jersey Bight subregions in the south-to-north direction in future decades, particularly in the midcentury. Future PSTs under WRF-HadGEM forcing are largely increased in the Hudson River as well as in west New England and east Long Island Sound counties for the 3-season PST. In the historical period, most WRF-HadGEM storms generating surge in the New England and Long Island Sound subregions passing to just offshore at southwest-to-northeast direction which are likely not as favorable to generating surge as storms that track further inland in the south-to-north direction. The WRF-CCSM4 RCM meanwhile shows more storms passing through or close to the subregions in the future decades than the historical one, and the storms are somewhat more intense. WRF-CCSM4 forcing showed increases to the 1-season PST in both future decades (greater increase in the midcentury) at the west New England and Hudson River counties. However, the WRF-CCSM4 storms are not as intense as some of those in future WRF-GFDL/WRF-HadGEM simulations, potentially explaining why the lower frequency 3-season PST under WRF-CCSM4 forcing is largely unchanged.
WRF-CCSM4 produces no particularly strong storm tide generating ETCs affecting the southern subregions (Delaware Bay, Chesapeake Bay, and North Carolina) in any decade (Figure 8). Regarding track patterns, there are generally more WRF-CCSM4 storm tracks passing through and just offshore of the North Carolina and the southern end of Chesapeake Bay in the historical decade than the future ones which could explain why future WRF-CCSM4 generated storm tides are moderately decreased (both the 3-season and 1-season PST) in these regions. WRF-GFDL produces a number of moderately intense storms directly impacting Delaware Bay and Chesapeake Bay in the historical decade. In future decades, WRF-GFDL storms are mostly less intense and there are fewer tracks that directly pass through or very close to the Delaware Bay and Chesapeake Bay subregions. However, WRF-GFDL storms tides are projected to be largely unchanged in these southern subregions with some small areas in Chesapeake Bay and North Carolina showing moderate increases. Noticeably, WRF-GFDL track angles in the historical decade are running very close to and parallel to the coastline while there is somewhat more variation in track angles for the future decades, which could be leading to certain tributaries of Chesapeake Bay having increased storm tides due to storm angles more favorable to surge. Furthermore, one may note the two moderately intense ETC tracks passing through the North Carolina subregion in the historical decade are in the west-to-east direction which would not be as favorable to surge in Albemarle and Pamlico Sounds as those tracks running just offshore of the subregion in the southwest-to-northeast direction plotted in the future decades. Compared to the historical decade, WRF-HadGEM storms affecting the Delaware Bay and Chesapeake Bay subregions are generally stronger and the track angle appears to be less oblique (more south-to-north running than parallel to the coast) in future decades. WRF-HadGEM generating storm tides in the northern areas Chesapeake Bay and Delaware Bay/River are indeed projected to increase in future decades while those in the areas exposed to the open ocean decrease. There also appear to be more intense WRF-HadG-EM storms tracking just offshore, as well as a couple just inland, of the North Carolina subregion in the future decades. Indeed, WRF-HadGEM generated storm tides are mostly larger in the North Carolina subregion in future decades.

Summary of Findings and Implications
Projected future changes in 1-season and 3-season peak storm tide elevations along the NEC under the RCP8.5 climate change scenario were found to range between ±0.3 and ±0.5 m, respectively. Variation due to RCM forcing is significant and generally greater than the variation between midcentury and late-century decades for the same RCM. In the New England and Long Island Sound subregions, there is no general consensus on midcentury or late-century changes to PST. This is similar to the findings of Lin et al. (2019) who project a small increase to cool-season surge heights under CCSM4 GCM forcing and a decrease under GFDL GCM forcing at Boston, MA (located in the New England subregion) for the mid-to-late 21st century . In particular, at Sussex County (#2) where Boston is located, our findings show little change to 1-season and 3-season PSTs under WRF-CCSM4 forcing and a larger decrease under WRF-GFDL forcing. Furthermore, our findings are largely in agreement to Roberts et al. (2017) and Lin et al. (2019) at New York County (#21) where we demonstrate little change to 1-season and 3-season PST under WRF-HadGEM and WRF-CCSM4 forcing and a decrease under WRF-GFDL forcing. Roberts et al. (2017) focused on New York City and found no significant change to the cool-season maximum surge elevation, while Lin et al. (2019) show that 1-3 season return period surge heights at New York City were slightly increased under CCSM4 GCM forcing and decreased under GFDL GCM forcing.
While there is not perfect consensus on projected changes to PST in the New York/New Jersey Bight, Delaware Bay, and Chesapeake Bay subregions, there is an indication that storm tides will increase at counties along the Hudson River, Delaware River, and northern end of Chesapeake Bay, and decrease at counties facing the open ocean along the mid-Atlantic Bight. This pattern is clearer for the late-century decade than the midcentury, and the magnitude of these changes varies fairly significantly between the RCMs. We found that increased PST in the upper reaches of the large estuaries and rivers could be attributed to more intense ETCs that track just inland (to the northeast) of these subregions in a south-to-north direction, where southerly winds to the southeast of these ETC centers will be favorable to generating surge in these locations. Zhang and Colle (2018) found that latent heating effects can cause ETCs to rapidly develop and become stronger in future warmer decades than historically, noting that this effect was more evident in high-resolution RCMs than coarse-resolution GCMs. Whereas, in the historical period, there appeared to be more ETCs tracking parallel to and just offshore of the coastline than in the future decades, likely resulting in the projected decreases to the PST for open coastal facing counties. To further that point, Figure 9 illustrates a unanimous reduction in ETC density over most of the ocean offshore of the NEC in the late-century decade, which is thought to be partly the result of weaker lower-tropospheric meridional temperature gradients in a warmer global climate (Chang, 2013;Colle et al., 2013). In contrast to the aforementioned storm tide projections, the 10-years surge height in northern Chesapeake Bay and Delaware Bay was projected to decrease for the 2054-2079 period while increasing along the mid-Atlantic Bight coast under CCSM4 and GFDL GCM forcing in Lin et al. (2019). This discrepancy could be at least partly due to the use of high-resolution RCMs in this study compared to the use of coarse-resolution GCMs, as well as other differences with the setups and analysis between the two studies (e.g., different time period of analysis, different return periods, storm surge instead of storm tide, and the storm surge heights were bias-corrected in Lin et al. (2019)). Furthermore, our finding for greater PSTs further inland was strongest under the WRF-HadGEM forcing which was not analyzed in Lin et al. (2019). There is little overall consensus on PST changes in the North Carolina subregion except for counties in the northern Albemarle Sound where the PST is projected to increase by all RCMs in both future decades. In these southern counties of the NEC region, the random tide-weather timing is not important, particularly in the Albemarle and Pamlico Sounds in which the tidal amplitude is smaller than 10 cm. Instead, storm tides are driven by passing ETCs locally generating surge in the sounds. In contrast, the larger tidal range in the New England, Long Island Sound, New York/New Jersey Bight, and Delaware Bay subregions leads to significant random uncertainty simply due to the phasing of the tides and weather conditions driving surge. Interestingly, although New England counties #1-5 located in the Gulf of Maine have the greatest tidal range in the NEC (see MHHW in Figure 5), the tide-weather timing related uncertainty is smaller than other northern counties such as those in Delaware Bay. This is likely related to the tide-surge interaction effects described in Horsburgh and Wilson (2007), which increase with the tidal range. It was found that modulation of the surge by the tide and vice-versa tends to result in the clustering of the peak nontidal residual and precludes it from coinciding with high tide, which in effect minimizes the impact of the random weather-tide timing on the peak storm tide level. Nevertheless, in most northern counties, including those in the Gulf of Maine, the direction (increase or decrease) of projected changes to 1-season and 3-season PST for each RCM forcing mostly depends on this random tide-weather timing. Alternatively, the magnitude of change could be much greater than the tidal phase-based average especially for the 3-season PST.
The importance of the aforementioned projected changes to the PST depend on the relative comparisons to the magnitude and uncertainty of future SLR. Assuming that the GCMs provide a reasonable uncertainty range of future SLR under the RCP8.5 scenario (see Section 4.2), projected SLR is 0.2-0.8 m by midcentury and 0.4-1.3 m by late-century. This implies that projected SLR and PST changes are about equally as important to consider under the low-end of SLR projections for midcentury. By late-century, low-end SLR will be slightly more important to consider for coastal flooding potential than any storm climatology-driven PST changes. Under the high-end SLR projection, even by midcentury potential PST changes are 2-4 times smaller in magnitude, and by late-century PST changes are 3-6 times smaller in magnitude. Although this study, as well as others Roberts et al., 2017), suggests that SLR will likely play a larger role in future changes to the cool-season coastal flooding potential in the 21st century, we should consider the combination of SLR and PSTs taking into account the full range of possibilities based on random tide-weather timing.

Uncertainties and Limitations
The WRF-based dynamically downscaled RCM simulations were conducted over fairly short time periods (decadal) and a relatively small number of GCM members (three) were used. For this reason, we avoided extrapolating our results to predict 100-years or other longer return periods using extreme value distributions. This study has instead focused on moderate frequency storm tide return periods (1 and 3 years) that ETC-driven surge events are mostly associated with in the NEC region. Furthermore, the RCM simulations were originally designed to investigate the North American continental climate and not particularly focused on resolving marine climatology. Future RCM simulations with greater computational resources will include larger portions of the ocean, more GCM members, and longer time periods.
On the hydrodynamic modeling side, the effects of wind-wave setup on coastal water elevations have been omitted in this study, primarily because wave modeling is significantly more computationally expensive than the hydrodynamic model. However, setup has been found to have a relatively small contribution to peak coastal water elevations in the NEC region (Marsooli & Lin, 2018), and is thus unlikely to impact our main finding. In addition, coastal flooding has been ignored, which if considered, generally results in lower water elevations on the ocean side compared to situations where inundation does not (or cannot) occur (Idier et al., 2019). Nevertheless, peak storm tide elevations recorded at the coast in our model simulations should generally be indicative of the coastal flooding potential.
We compared the magnitude of PST changes to SLR projections which we estimated from the parent GCMs for workflow self-consistency (avoiding external methodologies and models). It has been shown that CMIP5 GCMs may underestimate the externally driven anthropogenic component of SLR, particularly in the North Atlantic (Becker et al., 2016). Furthermore, widely varying representation of the Atlantic Meridional Overturning Circulation (AMOC) across GCMs is known to be a large source of uncertainty for SLR along the east coast of the United States (van de Wal et al., 2019). AMOC weakening has the potential to lead to rapid SLR for the NEC region (Yin et al., 2009). Compared to probabilistic SLR scenarios computed for the 21st century in the NEC region (Sweet et al., 2017), CCSM4/GFDL estimated SLR would indeed appear to correspond closest to the "low" scenario despite representing the RCP8.5 high-concentration pathway. However, HadGEM estimated SLR roughly corresponds to the "intermediate-high" scenario. We also note that storm tide dynamics and river discharge in the upper reaches of the estuaries and rivers may locally modulate SLR in a way that our analysis does not account for (Idier et al., 2019).

Data Availability Statement
Data sets for this research are available without restriction at Pringle (2020) under the Creative Commons Attribution 4.0 International license. National Oceanic and Atmospheric Administration (NOAA) tide gauge water elevation data were downloaded from https://tidesandcurrents.noaa.gov. National Buoy Data Center (NBDC) ocean buoy wind velocity data were downloaded from https://www.ndbc.noaa.gov. CCSM4 climate data were downloaded from https://www.earthsystemgrid.org/. GFDL-ESM2G and HadGEM2-ES climate data were downloaded from https://esgf-node.llnl.gov/search/cmip5. Specially, we chose the "historical" experiment for the historical decade and the "RCP8.5" experiment for the future decades, and the "r1i1p1" ensemble for both experiments. Geospatial data describing the boundaries of United States counties used in our analysis were downloaded from https://www.census.gov/geographies/mapping-files/time-series/geo/cartographic-boundary.html.