Understanding Tide Gauge Mean Sea Level Changes on the East Coast of North America

Sea level (SL) is an informative index of climate and a serious concern for coastal communities. Understanding the observational SL record is important from scientific and societal points of view. We consider the tide gauge SL record, focusing on data along the North American northeast coast, aiming to identify relevant geophysical processes responsible for observed SL changes. SL changes reflect dynamic and isostatic ocean effects. Recent works have interpreted accelerated and extreme SL changes along the northeast coast of North America primarily in terms of dynamic changes. In manuscript 1, we consider the influence of the ocean’s isostatic response to surface atmospheric pressure loading— the inverted barometer (IB) effect—on annual mean SL from tide gauge records. The IB effect explains ∼25% of interannual SL variance and accounts for ∼50% of the magnitude of a recent extreme event of SL rise along Atlantic Canada and New England. Estimated IB effects also amount to ∼10–30% of recent multidecadal SL accelerations over the Mid-Atlantic Bight and Southern New England. These findings reiterate the need for careful estimation and removal of isostatic effects for studies of dynamic SL. In manuscript 2, we continue our investigation of east coast tide gauge SL, seeking to better understand the relation between coastal SL and the variable ocean circulation. Annual SL records (adjusted for the IB effect) from tide gauges along the North American northeast coast over 1980–2010 are compared to a set of data-assimilating “ocean reanalysis” products as well as a global barotropic model solution forced with wind stress and barometric pressure. Correspondence between models and data depends strongly on model and location. At sites north of Cape Hatteras, the barotropic model shows as much (if not more) skill than ocean reanalyses, explaining ∼ 50% of the variance in the adjusted annual tide gauge SL records. Additional numerical experiments show that annual SL changes along this coast from the barotropic model are driven by local wind stress over the continental shelf and slope. This result is interpreted in the light of a simple dynamic framework, wherein bottom friction balances surface wind stress in the alongshore direction and geostrophy holds in the across-shore direction. Results highlight the importance of barotropic dynamics on coastal SL changes on interannual and decadal time scales; they also have implications for diagnosing errors in ocean reanalysis, using tide gauge records to infer past changes in ocean circulation, and identifying mechanisms responsible for projected regional SL rise. Finally, in manuscript 3, three global gridded reanalysis products are used alongside a carefully curated set of station records and tide gauges to consider IB changes over the global ocean and their relation to SL changes over the Twentieth Century. Centennial IB trends from reanalysis products show meridional structure consistent with the IB response expected under global warming. Annual IB variations show stronger amplitudes at higher latitudes, as in past studies focusing on higher frequencies or shorter periods. Discrepancies between gridded IB products tend to be smaller (larger) for more recent (earlier) periods and ocean regions with good (poor) historical data coverage. Comparisons between reanalysis products and station data reveal evidence for common errors across reanalyses over a wide range of time scales from annual to centennial. Notwithstanding their errors, the gridded reanalysis products are useful for interpreting tide gauge records: subtracting IB from SL records reduces both temporal variance within tide gauge records and spatial variance across tide gauge sites. Results advocate for making the IB correction to tide gauge SL data using reanalysis products in studies of ocean circulation and climate on centennial time scales.

gauge SL records. Additional numerical experiments show that annual SL changes along this coast from the barotropic model are driven by local wind stress over the continental shelf and slope. This result is interpreted in the light of a simple dynamic framework, wherein bottom friction balances surface wind stress in the alongshore direction and geostrophy holds in the across-shore direction. Results highlight the importance of barotropic dynamics on coastal SL changes on interannual and decadal time scales; they also have implications for diagnosing errors in ocean reanalysis, using tide gauge records to infer past changes in ocean circulation, and identifying mechanisms responsible for projected regional SL rise.
Finally, in manuscript 3, three global gridded reanalysis products are used alongside a carefully curated set of station records and tide gauges to consider IB changes over the global ocean and their relation to SL changes over the Twentieth Century. Centennial IB trends from reanalysis products show meridional structure consistent with the IB response expected under global warming. Annual IB variations show stronger amplitudes at higher latitudes, as in past studies focusing on higher frequencies or shorter periods. Discrepancies between gridded IB products tend to be smaller (larger) for more recent (earlier) periods and ocean regions with good (poor) historical data coverage. Comparisons between reanalysis products and station data reveal evidence for common errors across reanalyses over a wide range of time scales from annual to centennial. Notwithstanding their errors, the gridded reanalysis products are useful for interpreting tide gauge records: subtracting IB from SL records reduces both temporal variance within tide gauge records and spatial variance across tide gauge sites. Results advocate for making the IB correction to tide gauge SL data using reanalysis products in studies of ocean circulation and climate on centennial time scales.

ACKNOWLEDGMENTS
It's been a long road, and I've got lots of people to thank-Rui, meu professor, obrigado.
Kathy, advisor and tireless defender of my cause, I can't thank you enough.
Thank you also to my core committee members, Randy and Leonard.
To Annie, without whom I wouldn't have passed GFDII, thank you.
(I should probably thank NASA and the NSF for putting food on my table.) Thank you to my wonderful family, especially my darling sister and mother.
And, saving the best for last, to Hannah, my very best friend-Thank you. So much. This is yours as much as it's mine.  Table  Page 1.1 Tide gauges used here. Completeness is fraction of years for which data are available. Single daggers denote stations considered in the "northeast composite" of  that were used to produce the results in Figure 1.3. Double daggers denote stations considered in the "northeast hotspot" of  that were used to produce results in  1.3 (a) Interannual η from tide gauges (black) and η ib from different P a datasets (COREv2 blue; HadSLP2 red; NOAA 20CR yellow; NCEP-DOE R2 purple; ERA-Interim green) along New England andAtlantic Canada over 1979-2013 (in mm). Also indicated are magnitudes of the change over 2008-2010 (vertical lines) and η anomalies larger than the η standard deviation (dots). (b) Total interannual η variance explained by η ib (bars) and η variance explained by η ib when considering only the largest η anomalies as described in text (whiskers); (c) percent η ib contribution to the η change over 2008-2010; and (d)  Taylor diagram summarizing the correspondence between tide gauge records averaged north (circles) and south (squares) of Cape Hatteras and the corresponding sea level time series from the barotropic model (blue), GECCO2 (orange), ORAS4 (yellow), SODA (purple), and GODAS (green). Along the radial coordinate of the diagram is shown the standard deviation of the simulated ζ record divided by the standard deviation of the corresponding observational time series; along the azimuthal coordinate is shown the correlation coefficient r between the modeled and observed time series; and emanating from the reference point [i.e., the coordinate pair (1,1) starred on the diagram] is the relative root-mean-square deviation δ between the model and gauge records. The only significant correlation values are those from SODA and the barotropic model north of Cape Hatteras. (Note that the orange circle, corresponding to the performance of the GECCO2 product north of Cape Hatteras, is not missing from the figure, but rather falls outside the axis limits, owing to a negative correlation coefficient.  2.7 (a) Sea level from SHAL averaged over the 20 tide gauge sites north of Cape Hatteras (black) versus minus the average alongshore wind stress (here denoted τ ) over the northeast continental shelf (grey). Here we define the alongshore wind stress as the inner product between wind stress vector τ = (τ x , τ y ) and an alongshore unit vectorn = (cos ϑ, sin ϑ), where we have chosen ϑ = 30 • . We define the extent of the northeast continental shelf as the region between 53-100 • W and 35-45 • N where the ocean depth is < 1000 m. Sea level and alongshore wind on the northeast coast. Black curve in each panel is observed sea level record averaged over the 20 tide gauges north of Cape Hatteras (cf. Fig. 2.1a). Various colored curves in the different panels are the sea level predicted by averaging detrended annual alongshore wind stress anomalies over the shelf from various atmospheric reanalyses and scaling by −1 m 3 N −1 (see the text for more details): (a) NOAA-20CR (blue; , (b) ERA-20C (orange; , (c) ERA-Interim (yellow; Dee et al. 2011), (d) NCEP-NCAR (purple; Kalnay et al. 1996), and (e) NCEP-DOE (green; . We define alongshore wind stress and shelf extent as in Fig. 2 3.2 (a) Sample mean η ib linear trends T during 1900-2000 over the ocean across estimates (mm yr −1 ). Black lines are mean p a contours (10-hPa increments). Black squares and circles mark meteorological stations and tide gauges, respectively (cf. Fig. 3.1).

Introduction
Sea level changes can adversely impact coastal communities, leading to submergence, flooding, and erosion . Global mean sea level has risen steadily over the last century, with the rate of rise likely having accelerated during the last quarter-century . Moreover, magnitudes of extreme regional sea levels have likely increased over the last half-century [Rhein et al., 2013]. Given the flood exposure of growing coastal populations and assets , there is ample motivation to understand physical mechanisms underlying such past and present changes to project future sea levels and facilitate community adaptive planning [e.g., .
Regional sea level changes can reflect ocean dynamics, for example, adjustments of the gyre or overturning circulations to wind or buoyancy forcing; but they can also represent processes unrelated to ocean dynamics, such as the isostatic response to surface loading due to barometric pressure (or the inverted barometer effect) . Recent papers identify multidecadal accelerations Haigh et al., 2014;Sweet and Park, 2014; and extreme interannual anomalies  in tide gauge sea level along the northeast coast of North America; those works, and follow-on studies , interpret these regional sea level changes primarily in terms of ocean dynamics, for instance, changes in the Florida Current, Gulf Stream, meridional overturning circulation, as well as other coastal current regimes.
Given the focus of these papers, the impact of the inverted barometer effect on such multidecadal accelerations and interannual extremes in tide gauge records remains unclear. Some of these studies acknowledge and remove the inverted barometer prior to analysis, but without providing a sense of its impact, whereas a few others make no mention of this effect. Several of them briefly discuss the inverted barometer effect [e.g., Haigh et al., 2014;   , reiterating need for careful estimation and removal of isostatic effects for studies of dynamic sea level.

Datasets
To investigate sea level (η) between Virginia and Newfoundland, we use annual revised local reference time series from 30 tide gauges (  , we consider these records from 1950 onwards. Given our interest in multidecadal variations and interannual fluctuations, we use least squares to estimate and remove linear trends from the gauge records, which mainly reflect vertical land motion and global mean η changes . Consistent with past works , these tide gauge records evidence interannual and multidecadal η variations that are coherent across the Mid-Atlantic Bight, Gulf of Maine, and Atlantic Canada coast ( Figure 1.1). Recent studies attributing such η variations have highlighted the role of nearshore wind stress , baroclinic Rossby waves , overturning circulation changes , and divergence of Sverdrup transport .  20CR are 194820CR are -200920CR are , 185020CR are -201220CR are , 197920CR are -201420CR are , 197920CR are -201420CR are , and 187120CR are -2012tively. We compute spatial grids of η ib based on each of the different P a products according to where overbar is spatial average over the global ocean surface, g is acceleration due to gravity, and ρ is a constant reference ocean density. Annual η ib values are mapped to tide gauges using nearest neighbor interpolation only considering P a grid cells over the ocean. For most time periods, a η is not significant at the 1-σ level (Figure 1.2).

Multidecadal changes and accelerations
(We have adjusted the standard errors based on the autocorrelation of residuals [e.g., Haigh et al., 2014].) However, we find significant a η values over the most recent periods; for example, for 1950-2009, 1970-2009, and 1969-2011, we Figure 5). Averaging over periods with significant values, we find a mean a η value of 0.109 mm yr −2 .
We determine η ib contributions to the above η accelerations by computing the for all accessible combinations of P a dataset and time period with significant a η value. For a particular time period, we compute R based on a given P a product only if that product covers the entirety of that period. Contributions from η ib to a η are sensitive to choice of P a product and time period (Figure 1.2 et al. [2015] and this period is chosen as it is the longest one mostly covered by the various P a datasets. After removing quadratic fits to each time series over this period, we average tide gauge η records and the corresponding η ib estimates from the different P a products across the sites ( where we have defined the percentage of variance V in a time series x explained by another series y as (using σ 2 to denote variance) Considering the largest η anomalies, we notice that η ib contributes even more importantly; restricting focus to η anomalies with magnitudes larger than the sample standard deviation, the interannual η variance explained by η ib increases to 42% on average (Figure 1.3b).
The COREv2 η ib estimate explains noticeably less η variance than η ib from other P a datasets (

Discussion
We used tide gauge records and barometric pressure datasets to assess the This finding suggests that, while there are considerable uncertainties in available pressure datasets, the inverted barometer is not always negligible on multidecadal time scales, in broad agreement with . We also saw that the inverted barometer makes important interannual contributions, for example, amounting to to reconstruct past ocean current transports  and forecast future sea levels [Slangen et al., 2014].
This work complements recent studies of sea level at tide gauges along the northeast coast of North America that largely focus on dynamic changes rather than static effects. Although  and Haigh et al. [2014] regress Boston and New York sea level onto barometric pressure, winds are also included in their regression, so it is unclear how much the inverted barometer contributes to sea level at these cities.  notes that significant anticorrelation between the North Atlantic Oscillation and "nonlinear regional sea level anomaly" at Portland and Halifax is consistent with an inverted barometer, but he does not elaborate on how much this effect contributes at these sites. Our study builds on the foundation of these works, quantifying the inverted barometer's influence on accelerated multidecadal sea level rise along the northeast coast.
Woodworth et al. [2014] find that quadratic fits to tide gauge records along the North American northeast coast "are largely unaffected by applying the inverse barometer correction or not" over 1950and 1960. Similarly, Boon [2012 states that removing the inverted barometer from the Boston tide gauge "had minimal effect on the regression coefficients" of a quadratic fit to the data over 1969-2011. Our findings do not contradict those of  and , since one of our atmospheric reanalysis products also suggests a lesser role for the inverted barometer during these periods (Figure 1.2), but they do reveal that different pressure data can lead to distinct conclusions. Therefore, an important direction of future work, especially for interpretation of apparent sea level accelerations, will be better understanding of discrepancies between air pressure datasets at low frequencies, which could be due to differences in spatial resolution between products or other issues affecting reanalyses more generally or roughly three times larger than those found by , and indicate that isostatic effects were of central importance to this unique sea level rise event.
Our results underscore the need for careful consideration of the inverted barometer effect and any associated uncertainties for interpretation of tide gauge sea level records, lest the influence of isostatic effects be misinterpreted as a reflection of ocean dynamics.      that were used to produce the results in implications for diagnosing the uncertainties in current ocean reanalysis, using tide gauge records to infer past changes in ocean circulation, as well as identifying the physical mechanisms responsible for projected future regional sea level rise.

Introduction
Physical oceanographers have long sought to understand the relation between sea level on the northeast coast of North America and ocean dynamics in the North Atlantic. Appealing to simple models of the coastal response , earlier studies considered the connection between sea level fluctuations and local atmospheric forcing over the shallow continental shelf. Using two years of data,  reveals a strong link between adjusted sea level (that is, sea level corrected for the ocean's isostatic adjustment to barometric pressure changes) on the Nova Scotia shoreline and alongshore wind on the Scotian Shelf at periods > 20 days. This result is interpreted in light of a barotropic model, wherein the momentum balance is between wind stress and bottom drag in the alongshore direction, and geostrophic in the across-shore direction.  investigates sea level changes from long tide gauge records on the western boundary of the North Atlantic north of Cape Hatteras. The author hypothesizes that, while they are partly effected by local air pressure and wind stress, mean sea level anomalies along this coastline are also influenced changes in a wind-driven, coastally trapped boundary current.  contrast simulations from a homogeneous ocean model forced with air pressure and wind stress to tide gauge data on the North Atlantic western boundary. These authors discern that the model faithfully reproduces the observed adjusted sea level behavior on synoptic time scales (periods of 3-10 days).
This topic has also enjoyed renewed interest over the last decade, owing to concerns over global climate change, and the possibility that the ocean circulation will change and coastal sea level will rise (e.g., ). Based on geostrophic considerations and freshwater hosing experiments performed with a coarse resolution model,   While their findings are not necessarily contradictory, and may pertain strictly to particular time periods and frequency bands, the authors of these more recent dynamical studies are highlighting very different mechanisms in their interpretations of the tide gauge records. Yet for reconstructing past shifts in the ocean's general circulation  and anticipating future coastal sea level rise ), it is important to distinguish between the relative contributions of different ocean processes to sea level changes observed in the tide gauge record. With the goal of better understanding coastal sea level behavior, and partly motivated by Andres et al. (2013), who suggest the importance of barotropic dynamics, we study tide gauges, ocean reanalyses, and a barotropic model to address the following questions: • How well are year-to-year changes in sea level observed by tide gauges along the North American northeast coast reproduced by different ocean circulation models?
• Do barotropic processes contribute importantly to these observed sea level changes?
• What are the relative influences of local wind stress forcing over the shallow continental shelf versus remote wind driving over the deep open ocean?
The rest of this paper is structured as follows: in section 2, we describe methods and materials, namely the tide gauge data, ocean reanalysis products, and barotropic model solution; in section 3, we assess the skill of the ocean reanalyses and barotropic model in reproducing the tide gauge data; in section 4, model experiments are performed using the barotropic model to determine the roles of local and remote winds; finally, we conclude in section 5 with a discussion of our findings.

Methods and materials 2.3.1 Tide gauge records
To study sea level on the northeast coast of North America, we use annual revised local reference (RLR) records from 27 tide gauges (Table 1) We note that, while our main focus will be on tide gauges along the northeast coast, we have also included some tide gauges along the southeast coast of North America for purposes of comparison ( Fig. 2.1).
Here we focus on changes in dynamic sea level (ζ), hence we adjust the records for isostatic ocean response to barometric pressure (the inverted barometer effect), which can have an important impact on annual sea level changes in this area.
For example,  find that such air-pressure effects explain where overbar is spatial average over the ocean, g is gravity, and ρ is ocean density.
Values are mapped to gauge sites using nearest neighbor interpolation. Given our focus on ocean dynamics, we also remove estimated global mean sea level changes over the period .
Similar to recent works by  and , we restrict our focus to interannual and decadal changes. To isolate these time scales, we remove a linear trend from each of the annual tide gauge records.
This serves to filter out changes over longer periods due to global sea level rise and local vertical land motion ) and possibly also changes in thermohaline forcing and the Atlantic meridional overturning circulation   (1999).] In what follows, we seek to elucidate the dynamical mechanisms underlying these ζ fluctuations.

Ocean reanalysis products
To interpret the observed ζ anomalies (Fig. 2.1 Reanalyses were chosen largely based on their availability and temporal coverage. While each solution assimilates some ocean observations, two of them (GECCO2 and ORAS4) bring in altimetry data away from the coast, and none incorporate tide gauge data. A detailed description of the products is given below.
We take annual-mean ζ time series from the reanalyses. Since some models may not be faithful right at the coast, especially where the shelf is narrow compared to the model resolution, for each reanalysis and tide gauge, we map the model to the data by selecting the reanalysis ζ time series from the grid cell within a 300-

Barotropic model solution
To complement our study of tide gauge ζ based on ocean reanalyses, we also use a barotropic 3 model solution generated by the Massachusetts Institute of Technology general circulation model  The barotropic model setup makes use of a 900-s time stepping for the momentum equations along with a 3600-s time step for the free surface condition.
The model is started from rest using a 5-year spin-up period. During that time, it is driven with climatological P a and wind stress, whereafter it is forced with monthly reanalysis fields. While the model uses low-frequency (monthly) forcing, we also performed runs using high-frequency (daily) forcing fields, but they yielded nearly identical annual ζ solutions (not shown) and so are not discussed any further. To be consistent with the tide gauge records and ocean reanalyses, we remove the inverted barometer effect from the barotropic model solution. As with the reanalyses, we match model and data annual ζ fields by taking the nearby model ζ time series that explains the most variance in the tide gauge record. We remove a linear trend during the 1980-2010 period.

Comparing models and data
A number of recent papers compare tide gauges to sea level from ocean models in different areas ; . In order to gain deeper physical insight, we revisit this important topic, examining the tide gauge records and ocean model solutions along the North American northeast coast. To infer how well models reproduce the data, we compute two quantities: (1) the correlation coefficient (r) and (2) the relative root-mean-square deviation (δ) between the model and the data, where m and d represent model and data ζ time series, respectively, while σ is standard deviation.
The relationship between the models and the data varies from place to place and from model to model. There are no tide gauge sites at which the ζ data are significantly correlated with the modeled record from GECCO2, ORAS4, or GODAS ( Fig. 2.3b, c, e). Root-mean-square deviations between the data and either GECCO2 or GODAS are relatively large (greater than roughly 0.9; Fig. 2.3g, j).
ORAS4 performs only slightly better in this regard, for example, yielding δ ∼ 0.7 at Fernandina Beach ( Fig. 2.3h). These results are consistent with , who shows that GECCO2 has little skill in reproducing altimetric ζ data in this area In summary, our results show that ocean models differ in their ability to reproduce annual ζ changes observed on the North American east coast. They also suggest that barotropic processes contribute appreciably to interannual and decadal ζ variance on the coast north of Cape Hatteras. To elucidate the relevant barotropic dynamics, in the section that follows we report on results from additional numerical forcing simulations that were performed based on the barotropic model setup.

Forcing experiments and dynamical interpretation
Our simple barotropic model solution performs as well as, if not better than, other more complete (and data-assimilating) ocean general circulation model frameworks with regard to reproducing annual tide gauge observations along the northeast coast of North America. This demonstrates that more complex models do not necessarily produce more realistic solutions. In the most general terms, the ζ signals from the barotropic model can reflect dynamic ocean response to barometric pressure and wind stress locally as well as remotely. To reveal the roles of local and remote wind and pressure, we conduct the following experiments based on the barotropic model configuration: • PRES: In this experiment, we again run forward the barotropic model as described previously, but we "turn off" the wind stress surface forcing. Hence, once corrected for the inverted barometer effect, this solution represents the dynamic ocean response to barometric pressure.
• SHAL: For this run, we set to zero barometric pressure and wind stress over the deep ocean, leaving the wind stress over the shelf and slope (< 1000 m) as the only driver of ζ variability.
• DEEP: Similar to SHAL, we remove pressure and wind forcing over the shallow ocean from this simulation, allowing only wind stress over the deep ocean (> 1000 m) to force the model.
In all other respects (e.g., initial conditions), these perturbation runs are identical to the original barotropic ocean model simulation, which hereafter we refer to as the BASE experiment for clarity.
The outcomes of the experiments are summarized in Fig. 2 The PRES experiment evidences no appreciable dynamic behavior in this region and explains none of the ζ variance from the BASE simulation ( Fig. 2.6a). This result is not surprising, as the barotropic oceanic adjustment to pressure loading at these space and time scales is expected to be mainly isostatic and mostly explained by the inverted barometer response (e.g., . In sharp contrast, the ζ time series from the SHAL and BASE experiments are nearly identical-the correlation coefficient between them is 0.99 (Fig. 2.6b). This suggests that annual barotropic ζ fluctuations along the coast are driven by wind stress over the shelf and slope. The ζ fluctuations from the SHAL experiment are almost perfectly anticorrelated (coefficient of −0.99) with the local alongshore wind stress over the Mid-Atlantic Bight, Gulf of Maine, and Scotian Shelf (Fig. 2.7a).  also find strong anticorrelation between alongshore wind stress and coastal sea level, but the relation shown in Fig. 2.7a is much stronger than the one they see (cf. their Fig. 4b), likely because, as we use the barotropic component from the model rather than tide gauge data, we have effectively removed the influence of wind stress over the deep ocean and barometric pressure.  provides a physical framework for interpreting this antiphase relationship between sea level and alongshore wind stress. Consider a shelf of width W and depth H along the coast. Suppose that the momentum balance in the alongshore direction (hereŷ) is between wind stress and bottom friction, and say that geostrophy holds in the across-shore direction (x), whence, where f is the Coriolis parameter, g is the gravitational acceleration, v and τ y are the alongshore (i.e., meridional) velocity and wind stress, respectively, and A v is a vertical eddy viscosity. 4 If we assume that alongshore wind stress is constant, integrate across the shelf, and make substitutions with the equations, we obtain the following relation between sea level and alongshore wind stress,  we see in the DEEP simulation ( Fig. 2.7b), suggesting that barotropic Sverdrup balance is a plausible mechanism explaining this relationship.

Discussion
Previous investigations have studied the relation between coastal sea level and ocean circulation changes in observations of the past as well as projections of the future (e.g.,  although removing the barotropic model solution reduces the spectral power of the tide gauge data at all frequencies, the residual difference between them is slightly red. These findings accord with basic theory of the oceanic response (e.g., , which says that ocean stratification effects become more important with decreasing frequency. Additionally, the relationship between tide gauge and barotropic model sea level north of Cape Hatteras seems not to be stationary. For example, the correlation coefficient between these two time series is 0.91 for the decade 1983-1993 but 0.43 for the decade 1994-2004. Somewhat similarly,   to the annual tide gauge sea level records averaged over this coastline (Fig. 2.8).
We found that all alongshore wind stress products are significantly anticorrelated with the tide gauge records; after multiplying by the −1 m 3 N −1 scale factor determined in the last section, the reanalysis wind stress time series explain 44-55% of the annual variance in the tide gauge sea level record, depending on the choice of atmospheric reanalysis. This suggests that uncertainties in alongshore wind stress and local barotropic response are probably not responsible for the discrepancies in the skills of the different ocean models in this region (Fig. 2.5); rather, these discrepancies must be owing to inaccurate representation of some other forcing or process (e.g., thermohaline forcing, ocean stratification, baroclinic response, etc.).
Based on global analyses, Hernandez et al. (2014) and Balmaseda et al. (2015) find that models that assimilate altimetric data and have finer resolution generally reproduce tide gauge records better than solutions that either are more coarse or do not utilize altimetry. Thus, it might appear strange that the two models studied here that do incorporate altimetry (i.e., ORAS4 and GECCO2) perform poorly compared to other models that do not bring in this dataset (e.g., SODA). 5 However, it must be kept in mind (see below) that neither ORAS4 nor GECCO2 uses altimetric data near land. Notwithstanding concerns over potentially degraded quality of satellite altimetry data near the coast, the correspondence between standard altimetric products and tide gauge records can be good in some coastal regions (e.g., , and so it could be that the assimilation methods are discarding valuable data at the coast. Indeed, as specially tailored coastal altimetry products (e.g.,  come online and become more readily available, it will be important to bring them into ocean reanalyses for better representation of the coastal ocean. Another consideration is that representation of bathymetry could affect the along this stretch of coastline, which is small compared to the regional sea level rise anticipated during this coming century (e.g., Kopp et al. 2014;Slangen et al. 2014). We also found that, for a great majority (93%) of models considered, there is no significant change in the interannual alongshore wind stress variance over the duration of the simulation (not shown). Baltimore tide gauge records; these authors generally suggest that redistribution by oceanic Rossby or Kelvin waves may contribute to such regional sea level signals.
Analogously, based on a lagged correlation analysis considering European tide gauges,  suggest that westward wave propagation could result in multidecadal sea level oscillations at tide gauges between Halifax and Baltimore. However, it remains to determine how important variations in more local meteorological conditions are to multidecadal sea level changes along the coast. These important questions are beyond our current scope, and left for future study.

Description of ocean reanalysis products
The It employs the adjoint (or 4D-VAR) method to incorporate various satellite and in situ measurements including AVISO along-track ζ, mean dynamic topography, sea surface temperature from the AMSR-E satellite mission and the Hadley Centre Sea Ice and Sea Surface Temperature data set , as well as subsurface temperature and salinity from the EN3 database (Ingleby and Huddle ston 2007). Note that altimetric ζ fields are assimilated into the estimate only over regions deeper than 130 m. Bulk formulae are used for the adjusted surface forcing fields, which are based on the NCEP Reanalysis 1 (Kalnay et al. 1996;Kistler et al. 2001).

Acknowledgements
Author support came partly from NASA grant NNX14AJ51G. We thank Gaël  Table 1. Filled circles are correlation coefficients statistically significant at the 95% confidence level. Critical correlation coefficient values, determined for each pair of time series (von Storch and Zwiers 1999), are usually on the order 0.6-0.7. The black dashes separate sites north and south of Cape Hatteras.  Figure 2.7. (a) Sea level from SHAL averaged over the 20 tide gauge sites north of Cape Hatteras (black) versus minus the average alongshore wind stress (here denoted τ ) over the northeast continental shelf (grey). Here we define the alongshore wind stress as the inner product between wind stress vector τ = (τ x , τ y ) and an alongshore unit vectorn = (cos ϑ, sin ϑ), where we have chosen ϑ = 30 • . We define the extent of the northeast continental shelf as the region between 53-100 • W and 35-45 • N where the ocean depth is < 1000 m. (b) Sea level from DEEP averaged over the 20 tide gauge sites north of Cape Hatteras (black) versus wind stress curl integrated zonally across the deep ocean (> 1000 m) and averaged over 35-45 • N (grey). All the time series are detrended.  Kalnay et al. 1996), and (e) NCEP-DOE (green; . We define alongshore wind stress and shelf extent as in Fig. 2

Rhode Island
To be submitted to Journal of Atmospheric and Oceanic Technology

Abstract
Interpretation of tide gauge data in terms sea level (η) and ocean dynamics requires estimates of air pressure (p a ) to determine the ocean's isostatic responsethe inverted barometer effect (η ib ). Three gridded p a estimates (HadSLP2, NOAA-20CRv2, ERA-20C) are used alongside meteorological station p a and tide gauge η records to evaluate the contribution of η ib to η changes over the Twentieth Century. Agreement between gridded estimates is better during more recent periods and over regions with good historical data coverage, whereas it is worse for earlier time periods or in ocean areas with poor observational data coverage. Comparison against station data reveals the presence of systematic errors in the gridded estimates, for example, such that uncertainties estimated through differencing the gridded estimates underestimate the true errors by roughly 40% on interannual and decadal time scales. Notwithstanding such correlated errors, gridded estimates are still useful for interpretation of tide gauge data. Removing η ib estimates from η records reduces spatial variance in centennial trends across tide gauges by 10-30%, formal errors in centennial trends from individual gauges by ∼ 5%, and the temporal variance in detrended records by 10-15% on average (depending on choice of gridded estimate). Results here advocate for making the η ib correction to tide gauge records in studies of ocean circulation and global η over long, multidecadal and centennial time scales using an ensemble mean taken across several gridded η ib estimates.

Introduction
Tide gauges furnish some of the longest instrumental records of the ocean.
These records contain valuable signals related to global sea level and the ocean general circulation . However, tide gauge records must be interpreted carefully because they also measure the influences of processes unrelated to changes in the ocean's circulation or its volume . Surface loading owing to barometric pressure is one such process. Over periods longer than a few weeks, the oceanic adjustment to barometric pressure is more or less isostatic . The sea level responds like an "inverted barometer" in these cases , with water volume being redistributed in such a way that gradients in sea level cancel air pressure gradients, so that there is no signature in subsurface pressure, and no change in oceanic circulation.
To isolate the "signal" of changing ocean circulation or rising sea level, it is advisable to remove the "noise" related to the inverted barometer effect from tide gauge records. However, in numerous studies of global sea level  or the ocean general circulation  over the last century based on tide gauges, no correction is made for the inverted barometer effect. Some authors reason that this correction is small on century time scales ), while others argue that, due to increased scarcity of reliable air pressure data earlier in time, "no pressure adjustment is preferable to a partial or error ridden adjustment" .
Notwithstanding such concerns, the effects of air pressure on changes in sea level over very long (multidecadal and centennial) time periods, and the uncertainties associated with estimating them, have not been explored in any real systematic fashion. For example, some studies comment on the consistency between air pres-sure datasets or the impact of air pressure on sea level at a single tide gauge site ) or a small handful of locations ). Since they are based on a limited number of tide gauge records, or a single pressure dataset, it is difficult to generalize their findings. The quality of available pressure datasets and the relevance of the inverted barometer correction more generally as a function of space and time remains unclear.  studies the inverted barometer effect over the ocean during the period 1958-2000. He shows that pressure effects explain up to 40% of the variance in tide gauge records on monthly time scales, and that correcting for the inverted barometer reduces formal errors in sea level trends. Ponte recommends that, as more air pressure datasets become available, the impact of the inverted barometer correction on tide gauge records should be revisited for longer, centennial time periods.  find that air pressure effects made important contributions to recent sea level changes on the North American northeast coast , underscoring that such isostatic effects should be carefully removed in dynamical sea level studies. They also reveal that estimated inverted barometer contributions to sea level changes in this region on multidecadal time scales can be sensitive to the choice of air pressure dataset, which advocates for future investigations to better characterize uncertainties in these datasets over long time scales.
In light of the burgeoning number of pressure datasets covering the last century, it is timely to revisit the topic of sea level and the inverted barometer effect.
We address the following questions: 1. How consistent are air pressure reconstructions and atmospheric reanalysis products?
3. How much does the inverted barometer effect contribute to changes in sea level?
4. What is the best strategy for adjusting tide gauges for the inverted barometer effect?
The remainder of this study is structured as follows: in section 2, we describe the main datasets (i.e., meteorological station data, tide gauge records, and gridded reconstructions and atmospheric reanalyses); in section 3, we compare gridded air pressure estimates to evaluate their consistency; in section 4; we contrast those gridded estimates to a few rigorously vetted long meteorological station records to test for the presence of any systematic errors; in section 5, we evaluate inverted barometer effects on a carefully selected set of long tide gauge records; in section 6, we summarize the results, discussing the best practices for correcting tide gauge records for effects of air pressure.

Materials and methods
We investigate barometric pressure (p a ) changes and their contribution to sea level (η) changes. To this end, we use p a from reconstructions, atmospheric reanalyses, and meteorological stations, as well as η from tide gauges. We assess annual mean values over the common period 1900-2000.

Gridded p a products
We use three gridded European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis of the twentieth century (ERA-20C) from . These gridded estimates are described in more detail in section 3.8. In the broadest strokes, these products are generated by blending together gappy, sparse, discontinuous p a observational data based on advanced statistical methods. These products incorporate common datasets, and are not truly independent, but can differ in terms of analysis method, temporal coverage, horizontal resolution, quality control, bias correction, and detection of outliers. We evaluate the inverted barometer adjustment (η ib ) following Ponte (2006), where overbar is mean over the global ocean, ρ is ocean density, and g is gravitational acceleration.

Station data
We consider long p a time series from a few meteorological stations: Darwin, Chennai, Nagasaki, Azores, Gibraltar, Iceland, and Tahiti (Figure 3.1). These data have undergone intense scrutiny and quality control, which involves deleting dubious values, adjusting for step changes, comparing to paleoclimate data, and concatenating station records (e.g., Können et al. 2003). These records are taken from a NOAA-hosted website for the Working Group on Surface Pressure (http://www.esrl.noaa.gov/psd/gcos wgsp/Timeseries/) given with precision of 0.1 hPa.

Tide gauges
We also consider some long η records from coastal tide gauges ( Revised Local Reference (RLR) annual database (www.psmsl.org/data/). To avoid studying short, unrepresentative records, or ones that are overly influenced by solid earth processes, we applied stringent selection criteria to the tide gauges.
The 21

Changes in η ib from gridded products
In this section, we compare η ib changes over the global ocean during the Twentieth Century from the gridded reconstruction and reanalysis products to determine the consistency of these datasets. We find that discrepancies between gridded estimates, in terms of trends and variations, are smaller for more recent periods and over regions of the ocean having good historical data coverage, while they are larger for earlier time periods or areas with poorer historical data coverage ( Figures 3.2, 3.3).

Centennial trends
We compute η ib trends from the p a estimates using least squares. Let the η ib trend from estimate A ∈ {H, N, E} (where H, N , and E represent HadSLP2, NOAA-20CRv2, and ERA-20C) be T A . Suppose that T A represents the sum of the signal (true trend T S ) plus noise (error trend T A ), whence, The sample variance in trends across the estimates s 2 T then relates to the error trends according to, where n = 3 is the number of estimates and T . = n −1 A T A is the sample mean trend over estimates. In case of random errors, the sample variance χ a approaches the mean squared error trend χ b ; for common (positively correlated) errors across estimates, χ a represents a lower bound on χ b .
The sample mean trend T shows clear meridional structure (

Detrended fluctuations
Next, we examine the fluctuations in η ib that remain after removing the trends studied previously (Figure 3.2a-c). Similar to the trend analysis, we interpret each time series as signal plus noise, viz., can be related to error variances and covariances as, where brackets · are time mean. For uncorrelated errors, the left-hand side ψ a will tend toward the mean error variance ψ b . In the case of systematic errors, ψ a will be a lower bound on ψ b .
To with latitude has been seen in previous intraseasonal or interannual η ib or p a studies (e.g., . Such behavior is expected in case of quasi-geostrophic balance: given the change in Coriolis parameter with latitude, the magnitude of p a variation required to support a given wind fluctuation will increase towards the poles and decrease towards the equator. We evaluate √ ψ a over the ocean to see how well the estimates agree (Figure 3.2e). A latitudinal gradient is once again apparent. Differences between products are lower (< 5 mm) at low latitudes and higher (> 10 mm) at high latitudes. The smallest values (< 2 mm) are observed along the equatorial Indian Ocean, whereas the largest values (> 18 mm) are observed over the Pacific sector of the Southern Ocean.

point to the southeastern Pacific
Ocean as a region of especially pronounced uncertainty in the HadSLP2 product due to sparse historical data coverage.  shows similar patterns based on a comparison between annual fields from two reanalysis products over 1958-2000, highlighting the Pacific sector of the Southern Ocean as the area of largest differences between products. Compared to , we find comparatively larger differences between products (cf. our Figure 3.2e and his Figure 4b). This could be due to the fact that we compare more reanalysis products (cf. . Seeing as we consider a longer period extending more into the past, it could also reflect uncertainties earlier in the record or at lower frequencies.
Averaging over the global ocean suggests a typical error value of ∼ 5 mm.

Comparisons with station data
In this section, we compare the gridded reconstruction and reanalysis products to long records from meteorological stations to test for the presence of systematic errors in the gridded estimates. We find that errors can be correlated across gridded estimates and that some estimates perform better than others with respect to the data. In terms of trends, errors likely are correlated across the estimates, and HadSLP2 as likely as not performs better than either NOAA-20CRv2 or ERA-20C.
In terms of fluctuations, it is extremely likely that errors are correlated across the gridded estimates, and likely that HadSLP2 performs better than either NOAA-20CRv2 or ERA-20C (Figures 3.4-3.6).

Preliminaries
In the foregoing, it is difficult to partition errors between the estimates. Moreover, since products use similar data streams, they might share common errors.
A more thoroughgoing analysis would partition and separate errors in individual products (e.g., H 2 , N 2 , E 2 ) from covariance terms (e.g., H N , H E , N E ). One way to approach this problem (e.g.,  is to bring in a data source D with relatively small errors compared to the estimates, such that, e.g., (3.6) and, Our strategy here is to use data from a small number of meteorological stations that have undergone intense individual scrutiny (section 2.b), diagnosing all the terms in (3.3) and (3.5) at station locations. To see if these results are consistent with random normal errors (about a zero mean), we perform 100,000 trials of 21 draws (3 products × 7 sites) from the standard normal distribution (not shown). We interpret each of these random draws as the error from a particular product at a particular site. The probability that one product shows the smallest errors at four of the seven sites is P = 38%.
Similarly, we determine that one product showing smallest errors at more than three locations has a probability P = 52%. So, HadSLP2 showing smaller errors than both NOAA-20CRv2 and ERA-20C at more than half the stations is about as likely as not in the case of random errors. We also see that the probability that all three products show the same sign at three locations (more than two locations) is P = 17% (P = 24%). Therefore, we reason that the three products showing error trends of same sign at the Nagasaki, Azores, and Tahiti stations is unlikely given errors normally distributed about a zero mean. Finally, the probability of two products having errors of the same sign at all seven sites is P = 2.4%. Hence, we conclude that NOAA-20CRv2 and ERA-20C having error trends of the same sign at all stations is very unlikely given chance. Since these two estimates incorporate similar data (section 3.8), it is unsurprising that they show signs of common errors.

Detrended fluctuations
Are the error fluctuations [H , N , and E in (3.5)] correlated or uncorrelated across the estimates? To answer this question, we diagnose all the terms in the error budget (3.5) at the stations (Figure 3.6). We contrast those evaluations, based on comparisons between gridded products and station data, to an analogous exercise using 100,000 trials of 21 randomly simulated time series (not shown).
Each series has 101 entries ("1900-2000") randomly drawn from the standard normal distribution. Each time series is also scaled by a random factor drawn from the standard lognormal distribution. We interpret each of these series as representing the error fluctuations from a product at a station.
Strikingly, error fluctuations are positively correlated across the estimates at all sites (Figure 3.6). Such an occurrence is not seen in any of the 100,000 simulations of randomly generated series, and is therefore exceptionally unlikely in the case of random errors. At five of the seven stations, HadSLP2 shows the smallest error variance and NOAA-20CRv2 the largest error variance with respect to the station data. In our experiments with random series, the probability that one product has the smallest or largest variance at five (more than four) locations is P = 11% (P = 14% These findings demonstrate that error fluctuations covary across products. To Errors are also larger earlier in the record (Figure 3.4e-f).
The correlated errors strongly impact uncertainties estimated by comparing products (Figure 3.2e). For example, while there is a significant correlation (0.97) between ψ a and ψ b evaluated across the stations, the former underestimates the latter by a factor of ∼ 2 (Figure 3.6). This suggests that increasing the values in Figure 3.2e by ∼ 40% or so would paint a truer portrait of the average errors.
3.6 Impacts on tide gauge η records In this section, we compare the gridded reconstruction and reanalysis η ib products to η records at a few tide gauge sites to determine how much η ib contributes to η during the Twentieth Century. We find that removing estimates of η ib from η records reduces spatial variance in centennial trends across tide gauges by 10-30%, formal errors in centennial trends from individual gauge records by ∼ 5%, and temporal variance in detrended tide gauge records by 10-15% on average (Figures 3.7, 3.8).

Centennial trends
Centennial trends in η and η − η ib from gauge data and gridded products are shown in Figure 3.7a. On average, removing estimated η ib from η data reduces the variance in the linear trends across tide gauges. Reduction in spatial variance of the trend field ranges from ∼ 10% to 30% depending on choice of gridded product.
Moreover, formal error bars on trends in η − η ib generally tend to be somewhat smaller than formal errors in η trends. While formal trend error reduction can be > 10% at some locations (e.g., Brest, San Fransisco, Delfzijl, Newlyn), more generally it is ∼ 5%.
Given the limited number of records, these inferences could be very sensitive to this particular choice of tide gauges. To assess the sensitivity of results to tide gauge selection ("sampling error"), we perform 100,000 iterations of randomly choosing half of the records and recomputing values. In > 95% of cases, subtracting the ensemble mean η ib trend (averaged across the three estimates) from each gauge η record results in a reduction in the variance in the trends across the tide gauges.
Similarly, > 99% of the time, removing η ib from η reduces the average relative formal trend error. Thus, we conclude that results are robust and not sensitive to the particulars of tide gauge selection.

Detrended fluctuations
In Figure 3.7b, we show the percentage variance in each η record explained by the estimated η ib , where we define the percentage variance V in a time series x(t) explained by a time series y(t) as, where σ 2 denotes variance. The variance explained in η depends strongly on tide gauge location. Generally speaking, η ib tends to explain less η variance at lowerlatitude sites (e.g., Fernandina, Honolulu, Key West) and more variance at higherlatitude locations (e.g., Brest, Trieste, Newlyn). On average, η ib explains ∼ 10-15% of the η variance, depending on the choice of gridded product.
According to (3.8), η and η ib must agree in terms of amplitude and phase for η ib to explain high η variance, while low η variance explained by η ib result from differences in amplitude or phase.  (2015) show that η ib accounts for ∼ 25% of annual η variance during 1979-2013 on the North American northeast coast, we see that η ib explains little variance (< 10%) in annual η records along the east coast of the United States during 1900-2000.
Such contrasts between our results here and the findings of previous studies hint that the relation between η and η ib might not be stationary. Indeed, considering Newlyn ( Figure 3.8d), for example, it appears that the correspondence between η and η ib time series is stronger over 1950-1970 than during 1980-2000. To shine brighter light on the relationship between η and η ib as a function of time period and frequency band, we assess wavelet coherences ) at some sites ( with one another over the entire study period, whereas changes with periods ∼ 8 years are coherent more during the middle of the century. While wavelet coherence plots are similar for the two sites (cf. Figure 3.8e, h), η ib explains more η variance at Newlyn than at San Francisco since η ib amplitudes are more comparable to η magnitudes at the former than the latter place (Figure 3.8a, d).

Discussion
Tide gauges are one of the oldest observing systems of the ocean. For quantitative interpretation in terms of ocean dynamics or sea level, it is best to remove from tide gauge records the influence of isostatic volume redistribution related to surface pressure loading (the inverted barometer effect). We used gridded estimates and station data of surface air pressure (p a ) to investigate the inverted barometer effect (η ib ) and its impact on sea level (η) over the ocean during the Twentieth Century.
Centennial η ib trends (Figure 3.2a) evidence a meridional structure consistent with the observed strengthening of westerly circumpolar flows in the Southern Hemisphere that have been tied to stratospheric ozone depletion (e.g., . Detrended η ib fluctuations (Figure 3.2d) show stronger variability at higher latitudes and weaker variability at lower latitudes, as in previous intraseasonal and interannual studies , and as expected for quasi-geostrophic motions, given the change in the Coriolis parameter with latitude. The differences between the gridded η ib estimates are smaller for more recent periods and ocean regions with good historical data coverages, and are larger for earlier time periods or regions of the ocean with poor historical data coverage (Figures 3.2b, 3.2e, 3.3).
Comparisons against data records at meteorological stations show that errors can be correlated across gridded estimates and that some estimates perform better than others with respect to data. In terms of centennial trends, experiments with simulated random time series suggest that error trends likely are correlated across the gridded estimates considered here and also the HadSLP2 reconstruction as likely as not performs better with respect to the data than the NOAA-20CRv2 and ERA-20C reanalyses ( Figure 5). In terms of detrended fluctuations, it is extremely likely that error fluctuations are correlated across gridded estimates and likely that HadSLP2 performs better than NOAA-20CRv2 and ERA-20C ( Figure   6). Errors in NOAA-20CRv2 and ERA-20C are more strongly correlated with one another than with HadSLP2 ( Figure 3.4e-h, 3.5, 3.6), perhaps due to the fact that the former ingest such similar data streams, while the latter brings in many additional datasets not employed by the other products (section 3.8). These systematic errors can have a strong influence on uncertainties estimated through comparisons of gridded estimates. For example, in case of detrended fluctuations, variances in differences between gridded estimates underestimate the true error variances on those gridded estimates by roughly a factor of 2 on average ( Figure   6).
Notwithstanding these errors, gridded estimates are still useful for interpretation of tide gauges. Removing η ib estimates from η records reduces the spatial variance in centennial trends across tide gauges by 10-30%, the formal errors in centennial trends from individual tide gauges by ∼ 5%, and the temporal variance in the detrended tide gauge records by 10-15% on average (Figure 3.7). Therefore, while η ib contributions to η changes are not dominant on time scales considered here, consistent with past studies based on individual tide gauge records or single air pressure estimates (e.g., , neither are they completely negligible. Interestingly, the influence of η on η ib is a function not only of geographic location and frequency band, as noted in earlier works (e.g., , but also of time period (Figure 3.8). This finding elicits questions (beyond our scope) as to the reasons for the apparently non-stationary relationship between η and η ib .
Our results give guidance on how (and whether) to adjust tide gauge η records for the η ib effect. We advocate for making the η ib correction to tide gauge records in ocean circulation and global η studies on long, multidecadal and centennial time scales. While they can be characterized by (sometimes systematic) errors, gridded η ib estimates are of sufficiently good quality that adjusted tide gauge records typically have lower spatiotemporal variance that unadjusted tide gauge records.
Although uncertainties in η ib estimates can be large over some regions where historical p a data are scarce, long tide gauge η records tend to be found in regions where historical p a data are more plentiful, and so η ib estimates will tend to be relatively more well constrained at tide gauge sites. As a general "rule of thumb" for making the adjustment, we recommend using an ensemble mean across several gridded estimates. This approach, in addition to being objective, reduces variance within and across tide gauge records slightly more on average than individual estimates (Figure 3.7).  ]. This reconstruction is achieved based on reduced-space optimal interpolation (e.g.,    . Black lines are mean p a contours (10-hPa increments). Black squares and circles mark meteorological stations and tide gauges, respectively (cf. Fig. 3.1). (b) Square roots of trend sample variances √ χ a (mm yr −1 ). (c) Trend signal-to-noise ratios T / √ χ a (Fig. 3.2a/Fig. 3.2b). (d) Root mean η ib variance (RMV) over the ocean during 1900-2000 across products (mm).

Gridded p a products
(e) Values of √ ψ a (mm). (f ) Signal-to-noise ratios RMV/ √ ψ a for detrended fluctuations (Fig. 3.2d/Fig. 3.2e). Four white stars mark the locations of the η ib time series in Fig. 3.3. Thick dashed white contours denote a "cone of influence"  outside which edge effects may impact results.