Measuring global mean sea level changes with surface drifting buoys

Combining ocean model data and historical in-situ Lagrangian data, I show that an array of surface drifting buoys tracked by a Global Navigation Satellite System (GNSS), such as the GPS-tracked ones of the Global Drifter Program, could provide estimates of global mean sea level (GMSL) and its changes, including linear decadal trends. For a sustained array of 1250 globally distributed buoys with a standardized design, I demonstrate that GMSL decadal linear trend estimates with an uncertainty less than 0.3 mm yr$^{-1}$ could be achieved with GNSS daily random error of a couple meters in the vertical direction. This demonstration of the adequacy of the spatial sampling assumes that geodetic-grade GNSS vertical position measurements can be acquired from an array of drifting buoys, which is yet to be demonstrated. Development and implementation of such measurements could ultimately provide an independent and resilient observational system to infer natural and anthropogenic sea level changes, augmenting the on-going tide gauge and satellites records.


Introduction
Modern global mean sea level rise is an intrinsic measure of anthropogenic climate change. It is mainly the result of the thermal expansion of the warming ocean's water and the increase of ocean's mass from melting terrestrial ice [Rhein et al., 2013, Church et al., 2013. Global mean sea level rise is a major driver of the regional [Hamlington et al., 2018] and coastal sea level extremes Menéndez, 2015, Marcos andWoodworth, 2017] that impact millions of human lives and assets [Anthoff et al., 2006, Nicholls et al., 2011, and threatens ecosystems [Craft et al., 2009]. Currently, global mean and regional sea level changes are monitored by two observational systems: coastal and island tidal gauges [Douglas, 2001, Woodworth andPlayer, 2003] and satellite radar altimeters [Nerem, 1995, Ablain et al., 2017a. Tidal gauge records have high temporal resolution but their representativeness of the global mean sea level is biased towards the coasts. In contrast, the altimeter reference record is almost global (covering latitudes below 66 • ), but can only provide a near synoptic view about every 10 days. Here I show that if the satellite-tracked surface drifting buoys (hereafter drifters) of the Global Drifter Program recorded not only their geographical coordinates by 1 the Global Positioning System (GPS), but also their altitudes above the mean sea surface, the drifter array could provide estimates of global mean sea level changes, including longterm trends. With the current size of the drifter array, provided that the impact of potential biases is apprehended and quantified, global mean sea level decadal linear trend estimates with an uncertainty less than 0.3 mm yr −1 could be achieved with drifter GPS altitude random errors of a couple of meters. The drifter array could thus provide an independent and resilient observational system to infer natural and anthropogenic sea level changes, validating and augmenting the on-going tide gauge and satellites records.

Concept and Methods
In order to demonstrate the feasibility of the proposed idea, I take the example of NOAA's Global Drifter Program. Previously relying uniquely on the Argos system [Lumpkin and Pazos, 2007], this drifter array has almost completely transitioned to the Global Positioning System (GPS) for localization purpose, and to the Iridium satellite system to transmit their data to data centers on land [Elipot et al., 2016b]. Yet, despite the altitude coordinate being part of the positioning solution, only latitude and longitude data are retained in the transmitted data stream.

Concept
GPS altitude data at sea are a measure of height of the sea surface above a reference ellipsoid which center coincides with the Earth's center of mass. The relevant quantity is rather the height of the sea surface relative to a reference mean sea surface (MSS). The MSS is the sum of an equigeopotential surface called the geoid, which is slowly changing in time [Tamisiea, 2011] but is considered invariant for this study, and a mean dynamic topography associated with time-mean balanced oceanic motions over a reference time period. The interest here lies in relatively low-frequency changes of the mean sea level (MSL) which is the sea surface height (SSH) averaged over a given time period such as a day or a month. MSL is of particular interest when further averaged in space, regionally or globally, as its changes integrate physical processes associated with ocean dynamics and thermodynamics, such as its volume changes and its overall rise as a result of global climate change. As a consequence, GPS data acquired by a drifter would need to be subtracted from a co-located estimate of the reference MSS. The primary signal of interest is the local MSL where ε is an error, the sum of all physical contributions to SSH variability that are not relevant for global or regional MSL changes but may still affect their estimates. Such contributions include surface gravity waves, astronomically-forced tides, and internal gravity waves. The variance over time scales shorter than the time scale defining the MSL, in other words the uncertainty of h M SL that would be experienced by a drifter is since the MSS is constant in time and clearly independent of potential GPS heights from drifters. The purpose of this study is to simulate h M SL and quantify Var(h M SL ) for drifters, and subsequently calculate their global averages in order to assess if the drifter array could monitor GMSL changes with sufficient accuracy.

GPS errors
The main limitation of this study is that it is assumed that altitude measurements with sufficient accuracy can be acquired regularly from relatively small buoys drifting at sea, being ideally equipped with geodetic-grade GPS receivers (the GPS receivers currently in use have an overall estimated horizontal accuracy of 22 m [Elipot et al., 2016b]). This requires that enough computing power is present onboard each drifter to regularly calculate and transmit an accurate 3D position solution, or that time series of high frequency position data (typically at 1 Hz) could be transmitted intermittently to land for post-processing. Either solutions requires more electric power than currently available on a drifter to retain its 450-day designed life expectancy [Lumpkin et al., 2016]. In the rest of this study I assume that this power issue is solved, which is not unreasonable considering our current state of technological innovation. This study is therefore an ocean observing simulation, and is worthwhile as such, especially considering the critical importance of further measuring and understanding regional and global MSL changes. It is necessary to consider the systemic errors from GPS contributing to the terms Var(h GP S ) and 2Cov(h GP S , ε) in eq. (2). These errors include but are not limited to issues such as atmospheric effects, satellite orbit mismodeling, satellite and receiver clocks precision, and the choice of GPS calculation method [Santamaría-Gómez et al., 2011]. Despite such issues, GPS-equipped buoys at sea are successfully used to calibrate satellite radar altimeters, providing estimates of absolute biases of altimeter SSH data with an accuracy now reaching between 1 and 2 cm [Frappart et al., 2015, Born et al., 1994, Key et al., 1998, Watson et al., 2003. Historically, SSH measurements by GPS have relied on differential GPS techniques for which the GPS height of a buoy is calculated relative to a reference GPS station on land, typically no more than a few kilometers away [Watson et al., 2003, Watson et al., 2011. Such differential GPS method could not therefore be applied to the global drifter array for the purpose of estimating GMSL. Recently, new GPS methods relying on a single GPS receiver platforms have provided the opportunity to forgo reference land stations and therefore obtain SSH measurements in the absolute geocentric reference frame over the open ocean [Fund et al., 2013, Chen et al., 2013. In particular, the Precise Point Positioning (PPP) method can provide SSH measurements with accuracy as high as 5 cm when compared to 6-min tide gauge measurements [Kuo et al., 2012]. In fact, the PPP method has been successfully applied as an example to measure the surface height gradient of the Loch Ness with an unmanned water surface vehicle [Morales Maqueda et al., 2016], or to calibrate altimeter data with accuracy better than 2 cm [Frappart et al., 2015]. Whereas dedicated calibration experiments can rely on high frequency GPS data (typically at 1 Hz) and ample resources for post-deployment processing, the instrumental setup of drifters for geodetic-grade GPS measurements will likely be constrained by issues such as of battery ca-pacity, on-board processing power, and available bandwidth for transmitting position data. As a consequence, I do not attempt here to assign a unique value to the systemic individual GPS random error as a whole, but rather I consider a whole range of possible values, up to 2 m standard deviation for daily averages-a reasonable upper bound-and study the impact of all combined errors as a function of the number of drifters in the array.

Spatial sampling
Part of this study consists in verifying that the spatial and temporal sampling of the drifter array can provide GMSL estimates. For this I construct and analyze synthetic MSL measurements from hourly drifter trajectories [Elipot et al., 2016b] using sea surface height data from the GLORYS2V4 ocean reanalysis (Appendix A and Figure A1), excluding data in waters shallower than 120 m and poleward of 66 • N and S. This ocean model that notably assimilates altimetry data represents ocean dynamic sea level changes and global sea level rise [Ferry et al., 2012, Garric et al., 2017. Ocean models implicitly assume that the geoid and the reference ellipsoid coincide [Gregory et al., 2019] so that the time mean of the sea surface height of a model represents the mean dynamic topography portion only of the MSS. Therefore, using sea level anomaly (SLA) calculated as deviations from the time mean of the modeled sea surface height effectively simulates the term h GP S − h M SS in eq. (1). I linearly interpolate in time and space daily SLA from the GLORYS2V4 ocean reanalysis along drifter trajectories. I subsequently calculate a global mean at daily time steps to produce a relative time series of drifter GMSL changes (Appendix B). Does such a drifter GMSL time series represent true GMSL changes? To answer this question, this time series is compared next to a reference GMSL time series which is calculated from the ocean reanalysis using all grid points passing the same spatial selection criteria as the drifters.

Uncertainties from sampling error
The rise of GMSL induced by anthropogenic forcing is of primary importance, yet its estimate is still uncertain [Nerem et al., 2018, Kleinherenbrink et al., 2019. Can the drifter array capture the GMSL upward trend with sufficient accuracy despite its heterogeneous sampling? After subtracting seasonal cycle estimates (Appendix B), using daily time series, the linear trend of the drifter GMSL is 2.86 ± 0.03 mm yr −1 between 1993 and 2015.
Here, the uncertainty reported corresponds to 95% confidence interval calculated solely from the variance of the residuals from the least squares fit. For the same period, the reference GMSL trend is 3.30 mm yr −1 with an uncertainty less than 10 −2 mm yr −1 . The drifter GMSL therefore underestimates the true linear trends over 23 years by 0.44 mm yr −1 , but this is because of a large positive bias between 1993 and 1995 when the number of drifters is small (Figure 1a,c). Therefore, for the remainder of this paper I consider statistics and linear trends over the last decade of simulated data only, from the beginning of 2006 to the end of 2015. The year 2006 corresponds to the time when the drifter array reached maturity [Lumpkin et al., 2017] with approximately 1250 drifters ( Figure 1c). The drifter GMSL trend for 2006-2015 is 3.86 ± 0.03 mm yr −1 , not significantly different from the reference trend of 3.82 ± 0.02 mm yr −1 . Both estimates show the most recent acceleration of GMSL rise [Nerem et al., 2018] rendered by the ocean reanalysis. In conclusion, the absolute difference of decadal trends, and thus the expected inaccuracy of GMSL trend because of sampling errors only is 0.04 mm yr −1 . For the period 2006-2015, the rms difference between the GMSL daily time series, after subtracting seasonal cycle estimates, is 4.4 mm. This error for GMSL, not formally represented in eq. (2) because it arises from the global averaging calculation, is effectively the sampling error of the drifter array, which will be ultimately counted independently towards the total error. Such GMSL error implies a GMSL trend error of 0.03 mm yr −1 and is therefore nearly consistent with the trend error just diagnosed directly.
The Global Climate Observing System (GCOS) recommends in its 2016 implementation plan [Belward, 2016] that GMSL and GMSL decadal trend, two Essential Climate Variables, be reported on weekly to monthly time scales with an absolute stability, or uncertainty, of 2-4 mm and 0.3 mm yr −1 , respectively. I investigate next if such levels of accuracy can be achieved when adding to the drifter GMSL sampling error, just considered above, other possible sources of errors, systematic or random in character.

Error budget and analysis 4.1 Random errors
Variance in daily MSL from hourly heights acquired by drifters will arise from relevant physical processes not represented by the model of the GLORYS2V4 ocean reanalysis: SSH variability from barotropic tides, from baroclinic tides, and from surface gravity waves.
SSH variability from barotropic tides and baroclinic tides are estimated here using maps of SSH variance from a global 1-year run of the HYCOM numerical model forced by atmospheric fields and tidal potential [Savage et al., 2017]. Specifically, SSH variances from the diurnal, semi-diurnal and supertidal frequency bands are interpolated at all drifter hourly locations ( Figure A2, median value of corresponding standard deviations is 0.19 m), and subsequently globally averaged to derive a daily time series of error ( Figure 1b). This error decreases steadily as the number of drifters increases, and becomes relatively stable after 2006, averaging to 1.6 mm between 2006 and 2015. When compared to a discrete number of in-situ observations, HYCOM can underestimate the observed high-frequency SSH variance [Savage et al., 2017], thus the MSL errors from these processes may be underestimated here. However, these estimates may still be reasonable since a fraction of the stationary component of the SSH variability in tidal bands could be subtracted from drifter observations using predictions from a barotropic tidal model [Egbert and Erofeeva, 2002], and the error used here could be actually reduced. In addition, on-going research on the stationary component of baroclinic waves suggest that part of their SSH signal may also be subtracted [Zaron andRay, 2017, Zaron, 2019].
SSH variability from surface gravity waves is estimated using significant wave height (SWH) data obtained from a global, 0.5 • spatial resolution, 3-hr temporal resolution, run of the wave model WAVEWATCH III as part of the IOWAGA project (Appendix A), forced by National Centers for Environmental Prediction Climate Forecast System Reanalysis winds [Ardhuin et al., 2011]. The SSH variance fields from surface gravity waves calculated as (SWH) 2 /16 are linearly interpolated at the drifter locations at hourly time steps along their trajectories ( Figure A2, median value of corresponding standard deviations is 0.58 m). When globally averaged to form a daily time series, the GMSL error from surface gravity waves also decreases as the number of drifters increases and becomes stable after 2006, averaging to 4.4 mm for 2006-2015 (Figure 2b). The next SSH variance error arises from uncertainty in the MSS subtraction (eqs. 1 and 2), and is quantified using the formal error of the CNES CLS15 MSS product at 1-minute spatial resolution [Schaeffer et al., 2012]. The distribution of MSS errors sampled by the drifters ranges from 2 mm to 0.42 m ( Figure A1), but averages to 0.1 mm for the global mean for 2006-2015 (Figure 1b).
The three SSH error time series just derived above all anti-correlate strongly with the number of drifters time series (Figure 1c), as well as does their sum. This anti-correlation makes it possible ultimately to predict how the total random error can be expected to further decrease linearly if the number of drifters is increased.

Bias errors
Arguably, the formal MSS error considered earlier may be an underestimate of the total MSS error that also includes modeling or methodological biases [Pujol et al., 2018]. In order to quantify this additional error, I compare two MSS products, CNES CLS15 and DTU15 [Andersen et al., 2015], which use overlapping altimetry datasets but use different methodologies. I linearly interpolate the time-invariant but spatially-varying difference of MSS (CLS15 minus DTU15) along drifter trajectories ( Figure A2, 8.65 mm median standard deviation), and calculate a global average at daily time steps to produce a time series of MSS bias (Figure 3a). This bias time series effectively induces an error that takes the form of added variance for the drifter GMSL estimate. The corresponding error amounts to 0.55 mm on average for 2006-2015.
The next systematic source of errors to consider is an instrument-to-instrument bias due to the positioning of the GPS antenna within the buoy of an individual drifter. To measure sea level, the use of a buoyant float containing a GPS antenna requires to evaluate the vertical distance between the flotation line and the GPS antenna phase center. Using reference tide gauge data, calibration of GPS antenna on buoys can provide such distance with an accuracy better than 1 cm [Frappart et al., 2015]. However, in the case of the drifter array, calibration of all buoys would not be practically feasible so that an uncertainty will arise from unavoidable deviations in the construction of each drifter, implying a constant bias for the life of a single drifter. I assess the potential impact of such instrument-toinstrument bias by assigning to each individual drifter a random constant bias drawn from a normal distribution with a standard deviation of 1 cm, and subsequently recalculate the drifter GMSL. This operation is conducted one hundred times and the variance of all these 6 alternate GMSL time series is calculated at daily time steps (Figure 2c). This variance effectively takes the form of an added random error, averaging 0.32 mm for 2006-2015. This variance is however anti-correlated (at -0.67) with the number of drifters which implies that it can be expected to reduce with an increased number of drifters.
Another possible source of errors arises from the fact that the surface float of a drogued drifter has been observed to sink underwater as wave crests passe its location, being pulled underwater by the tethered drogue [R. Lumpkin, personal communication]. As the height of a drifter is eventually calculated by the GPS from data acquired over a finite time window every hour (the exact processing varies between the different drifter manufacturers), the sinking behavior of the buoy could introduce a negative bias as a drogued drifter would sample preferentially the troughs over the crests of surface gravity waves. I quantify this potential bias by considering again the outputs of a numerical simulation of waves [Ardhuin et al., 2011], but this time taking the standard deviations of the wave field as being the scale parameters for probability distribution functions which are originally normal and centered-a reasonable first approximation for SLA induced by waves [Longuet-Higgins, 1963]-but having no value above one standard deviation. The new means of such truncated distributions effectively provide the biases which are linearly interpolated in time and space at drogued drifter positions only. The wave bias is ignored for undrogued drifters as these ones are expected to ride the waves. The true in situ behaviors of drifters are actually unknown, but the two extreme behaviors considered here should be representative of the range of drifter behaviors. The impact of these biases on trend estimates suggest that in the presence of wave biases, it will be necessary to calculate separately GMSL estimates for drogued and undrogued drifters. Indeed, averaging simulated data from both biased drogued drifters and unbiased undrogued drifters leads to an overall GMSL fall, rather than a rise, during for the 2006-2015 period (Figure 2b). This occurs because during approximately the years 2013 and 2014 the number of drogued drifters surpassed the number of undrogued drifters (Figure 1c). In contrast, using unbiased undrogued drifters only (on average 670 daily over the last decade) provides a GMSL decadal linear trend of 3.72 ± 0.03 mm yr −1 , an error of only -0.10 mm yr −1 compared to the reference trend. Using biased drogued drifters only (on average 394 daily) provides a GMSL decadal linear trend of 4.26 ± 0.07 mm yr −1 , an error of 0.44 mm yr −1 . Ultimately, the impact the wave bias becomes an issue of sampling bias and the number of drifters to be used to reliably estimate trends, considered below.

Error budget
As a final step, I combine all GMSL errors discussed previously to derive estimates of daily and monthly GMSL uncertainty and decadal trend uncertainty based on daily values, as function of the number of drifters and a range of GPS random error. The error variance of the GMSL daily estimates for the 2006-2015 time period is estimated as the sum of five terms: Next, regression theory predicts that the error variance of GMSL propagates into the error variance for the linear trend parameter β 1 as where the t k are the daily time steps for the time period over which β 1 is estimated and t is the mid-point in time.
The sampling error variance, σ 2 S , is a constant equal to the rms value of the difference between the drifter and reference GMSL estimates for 2006-2015, after subtracting seasonal cycles (σ S = 4.40 mm). σ 2 MSS bias is the variance from the MSS bias (Figure 2a), taken as a constant equal to the average for 2006-2015 (σ MSS = 0.55 mm). σ 2 antenna is the variance from the antenna bias which is estimated at every daily time steps (Figure 2c), averaging to 0.32 mm for 2006-2015. This last variance is anti-correlated with the number of drifters N (correlation -0.67). The SSH variance from unresolved ocean physics, σ 2 SSH , estimated at daily time steps is also strongly anti-correlated with N (correlation -0.90). Thus, I derive the linear model σ 2 antenna + σ 2 SSH = α 1 N + α 2 with α 1 = −0.0203 mm 2 per drifter and α 2 = 6.6 2 mm 2 . In order for such model not to predict negative variance, it is set to zero for all values of N larger than 2160. The last term of eq. (3) is a simplification of a global average calculation (Appendix B) that neglects the area-weighting. This term also assumes that the GPS height daily-averaged error variance σ 2 GPS is the same for all drifter hourly measurements. Further, when used to propagate into trend errors, it is assumed that these globally averaged GPS errors are uncorrelated day to day. GPS errors from land stations are typically correlated [Mao et al., 1999, Zhang et al., 1997, their spectrum being often characterized as a mixture of low-frequency white noise and high-frequency power law noise [Santamaría-Gómez et al., 2011]. Nevertheless, I assume that daily averaging of hourly measurements of a single drifter, as well as globally averaging over all drifters, should reasonably filter out correlated errors for GMSL. This assumption will only be testable once surface drifters record and transmit their altitude data.
Yet, as an exercise, varying the values of N and σ GPS produces a range of estimates for σ GMSL (Figure 3a). The same calculation is re-conducted for monthly GMSL values (Appendix B and Figure 3b). Finally, the uncertainty for GMSL decadal linear trend based on daily errors is calculated using eq. (4) (Figure 3c). For the nominal size of the drifter array (1250 drifters), however small the magnitude of the GPS error is, the uncertainty for GMSL daily estimates is always larger than 6 mm and increases rapidly with the GPS error. For monthly GMSL estimates, if the standard deviation of daily GPS error is less than 40 cm, the uncertainty remains below 4 mm which is the upper bound of the GMSL errors recommended by GCOS (2 to 4 mm). As for an estimate of the GMSL decadal linear trend based on 1250 drifters, this one would be less than 0.3 mm yr −1 , the upper value recommended by GCOS, if the GPS daily error remains below 1.8 m standard deviation. If only undrogued drifters unaffected by a potential wave bias are used (670 on average for the [2006][2007][2008][2009][2010][2011][2012][2013][2014][2015], the uncertainty for the decadal linear trend would be below 0.3 mm yr −1 if the GPS daily error is maintained below 1.3 m. These errors are more than an order of magnitude larger than the reported errors for geodetic buoys at sea using PPP GPS techniques [Fund et al., 2013, Frappart et al., 2015. As a point of comparison, It is estimated that the uncertainty for GMSL long-term trend estimates from altimetry remains on the order of 0.5 mm yr −1 because of atmospheric signal corrections, satellite orbit uncertainties, altimeter instrument drifts, and processing methodologies [Ablain et al., 2009, Ablain et al., 2017b.

Conclusions
In conclusion, this study aimed to demonstrate that the global drifter array could provide a third means-in addition to satellite altimetry and tide gauges-to monitor GMSL changes related to climate processes, provided that drifters are tracked using geodetic-grade GNSS technologies. A third observational system to measure mean sea level would be greatly beneficial to further validate and calibrate the existing observing systems, as well as refine the closing of the sea-level budget and advance knowledge. Importantly, drifters could provide a reliable estimate of GMSL trends related to climate change. An estimate of GMSL from drifters would not have the geographical limitation of the tide gauge network, and might provide truly synoptic GMSL estimates at higher frequency than derived from the reference altimetric satellites. The ocean observing system simulated in this study is currently not achievable because of the anticipated limitation of the GPS receivers currently equipping the drifters of the global array. This limitation cannot be fully tested because the position records of GPS-tracked surface drifters do not currently contain altitude. Considering the median lifetime of 240 days of drifters between 2006 and 2015, if the Global Drifter Program array changed its specification to equip drifters with standardized GPS receivers with a tested daily random error of a few meters in the vertical direction, and provided that potential biases are apprehended and quantified, the drifter array could be providing its first GMSL estimates within a few years, augmenting the on-going satellite and tide gauge observational networks. 1993 1995 1997 1999 2001 2003 2005 2007 2009 Figure 1: Global mean sea level estimates. a. GMSL time series estimates. The blue line is the estimate from simulated drifter data and the red line is the reference estimate from the GLORYS2V4 ocean reanalysis. Both estimates exclude data in water depth less than 120 m and poleward of 66 • north and south. The thicker lines are the 9-month smooth estimates after removal of seasonal cycle estimates. b. Daily error estimates for the drifter GMSL time series. The sampling error is a constant equal to the RMS of differences between the deseasoned daily GMSL times series from GLORYS2V4 and the drifter array. c. Number of drifters per day in the hourly drifter database.

Appendices A Data sources
The hourly drifter database [Elipot et al., 2016a] is available from the website of the Data Assembly Center of the Global Drifter Program, hosted at the NOAA/Atlantic Oceanographic and Meteorological Laboratory (http://www.aoml.noaa.gov/phod/dac/hourly data.php). Out of the database version 1.01, I selected 12,784 unique drifters between January 1, 1993 and December 31, 2015, totaling 115,344,433 hourly locations of drifters both drogued and undrogued. The hourly database contains only continuous drifter trajectory segments of 12-hr length or longer. Sea surface height data used for simulating drifter sea surface height measurements are from GLORYS2V4, an ocean reanalysis product for the reference altimetry era, 1993 to present (Figure ??). The GLORYS2V4 reanalysis is based on an ocean and sea-ice general circulation model (NEMO) with 1/4 • horizontal resolution, and assimilates sea surface temperature, in situ profiles of temperature and salinity, and along-track sea level anomaly observations from multiple satellite altimeters [Ferry et al., 2012]. GLORYS2V4 captures climate signals and trends, and realistically represents mesoscale variability [Garric et al., 2017]. For this study, the global-analysis-forecast-phy-001-025 files of daily means of sea surface height were downloaded after registration from the Copernicus Marine Environment Monitoring Service (CMEMS) at http://marine.copernicus.eu. Sea level anomalies are calculated by subtracting the time-average of the model sea surface heights over the 1993-2015 time period.
Time-mean maps of sea surface height variance from barotropic tides and internal gravity waves at diurnal, semi-diurnal and supertidal frequencies from the HYCOM numerical model were provided as electronic files by Anna Savage [Savage et al., 2017]. The maps available on the original tri-pole grid of the HYCOM model were recasted in approximate 1/4 • bins before spatial interpolation onto the drifter locations.
Significant wave height (SWH) data are from the wave model WAVEWATCH III from a global run at 0.5 • spatial resolution and 3-hr temporal resolution, forced by National Centers for Environmental Prediction Climate Forecast System Reanalysis winds (NCEP CFSR), as part of the IOWAGA project (https://wwz.ifremer.fr/iowaga/, [Ardhuin et al., 2011]). These data for 1993-2015 were downloaded from the IOWAGA ftp site (ftp://ftp.ifremer.fr/ifremer/ww3/HINDCAST/). The SSH variance fields from surface gravity waves calculated as σ 2 = (SWH) 2 /16 were subsequently interpolated linearly in time and space at the drifter locations at hourly time steps along their trajectories.
The CNES CLS15 MSS product at 1-minute resolution including and its formal error field are distributed by AVISO (https://www.aviso.altimetry.fr/en/data/products/auxiliaryproducts/mss.html). This MSS product has been derived using 20 years of data from 1993 to 2012, as an update of the CNES CLS11 product [Schaeffer et al., 2012]. The formal error from the optimal interpolation mapping technique takes into account altimetric noise, ocean variability noise, and along-track biases.
The DTU15 MSS and its formal error at 1-minute resolution [Andersen et al., 2015] were obtained from the ftp server of the Danish National Space Center (ftp://ftp.space.dtu.dk).

B Details of methods
A global or regional mean value x, from individual measurements x i of a variable x is calculated as an area-weighted mean as where w i = cos φ i with φ i the latitude of measurement x i [Henry et al., 2014, Masters et al., 2012. This formula is used to compute GMSL at daily time steps from simulated sea surface height data from drifters, selecting data equatorward of 66 • and in waters deeper than 120 m. Another method to compute MSL consists in first averaging data in regular geographical bins, and second averaging all bin values with area-weighting (with the weights being the cosine of the latitude of bin centers) [Henry et al., 2014]. This second method has been tested and the results are qualitatively the same. If σ 2 i is the independent error variance of individual measurement x i , then the variance of x defined by equation (5) is the square root of which provides the standard error of x. The variances of the drifter GMSL at daily time steps are calculated using this formula with σ 2 i taken as the estimated individual variances from unresolved physical processed (surface gravity waves, tides, and internal waves) and the formal CLS15 MSS error variances.
Linear trend estimates and their corresponding standard error estimates for SL times series are calculated by ordinary least squares for which residuals from the fit are used as estimate for the variance of the data.
In order to be eventually removed, the seasonal cycle of the GMSL time series are estimated by bandpass filtering. The time series are convolved with a Hanning window of 2-year length multiplied by complex exponentials at the annual and triannual frequencies. Mirror boundary conditions are applied at the end of the time series to limit the impact of edge effects. The frequencies for the seasonal filters are chosen based on the 95% threshold of false alarm for the coherence squared between the drifter-based GMSL time series and the model GMSL time series (not shown). The smoothing filter used to extract monthly and interannual variability of the GMSL time series is a NadarayaWatson kernel estimator with an Epanechnikov kernel of half-width 1 month and 9 months, respectively [Fan and Gijbels, 1996]. Figure A1: Mean Sea level characteristics of the GLORYS2V4 reanalysis and density of hourly drifter observations. a. Linear sea level trends from 1993 to 2015 of GLORYS2V4 minus the GMSL trend of 3.30 mm yr −1 . b. Seasonal amplitude of mean sea level from GLORYS2V4. c. Drifter hourly observations density in 1/4 • bins on average for each year since 2006. Land masses are shaded white. Figure A2: Distribution of error estimates. Histograms of interpolated mean sea level errors at all drifter locations (square root of interpolated variances) from surface gravity waves, from the formal errors of the CNES CLS15 mean sea surface, and from tides and internal waves. The histogram of absolute differences between the CNES CLS15 and DTU15 mean sea surface products is also shown. Vertical dashed lines indicate the median value of each populations.