Sub-basin scale sea level budgets from satellite altimetry, Argo floats and satellite gravimetry in the North Atlantic

Abstract. In this study for the first time an attempt is made to close the sea level budget on a sub-basin scale in terms of trend, annual amplitude and residual time series, after removing the trend, the semi-annual and annual signals. To obtain errors for altimetry and Argo full variance-covariance matrices are computed using correlation functions and their errors are fully propagated. For altimetry we apply a geographically dependent intermission bias (Ablain et al., 2015), which leads to differences in trends up to 0.8 mm yr−1. Since Argo float measurements are non-homogeneously spaced, steric sea levels 5 are first objectively interpolated onto a grid before averaging. For the Gravity Recovery And Climate Experiment (GRACE) gravity fields full variance-covariance matrices are used to propagate errors and statistically filter the gravity fields. We use four different filtered gravity field solutions and determine which post-processing strategy is best for budget closure. As a reference the standard 96-degree DDK5-filtered CSR solution is used to compute OBP. A comparison is made with two anistropic Wiener-filtered CSR solutions up to d/o 60 and 96 and a Wiener-filtered 90-degree ITSG solution. Budgets are computed for 10 ten polygons in the North Atlantic, defined in a way that the error on the trend of Ocean Bottom Pressure (OBP) + steric sea level remains within 1 mm yr−1. Using the anisotropic Wiener filter on CSR gravity fields expanded up to spherical harmonic degree 96, it is possible to close the sea level budget in nine-out-of-ten sub-basins in terms of trend. Wiener-filtered ITSG and the standard DDK5-filtered CSR solutions also close the trend budget, if a Glacial Isostatic Adjustment (GIA) correction error of 10-20 % is applied, however the performance of the DDK5-filtered solution strongly depends on the orientation of the polygon 15 due to residual striping. In seven-out-of-ten sub-basins the budget of the annual cycle is closed, using the DDK5-filtered CSR or the Wiener-filtered ITSG solutions. The Wiener-filtered 60and 96-degree CSR solution in combination with Argo lack amplitude and suffer from what appears to be hydrological leakage in the Amazon and Sahel regions. After reducing the trend, semi-annual and annual signals, 24-53 % of the residual variance in altimetry-derived sea level time series is explained by the combination of Argo steric sea levels and Wiener-filtered ITSG OBP. Based on this, we believe that the best overall solution 20 for the OBP component of the sub-basin scale budgets is the Wiener-filtered ITSG gravity fields. The interannual variability is primarily a steric signal in the North Atlantic, so for this the choice of filter and gravity field solution is not really significant.


Introduction
Several studies attempted to close the sea level budget by using satellite altimetry, satellite gravimetry and observations or reanalyses of ocean temperature and salinity on a global scale. Closure of the budgets is required to get a consistent division between mass or Ocean Bottom Pressure (OBP) and steric related sea level changes. This helps us to identify the contributors to present day sea level changes. Contributors that affect ocean OBP are glacier and ice sheet melt and land water storage, 5 while heat fluxes between ocean and atmosphere contribute to steric changes. Ocean dynamics have an effect on both the OBP and the steric change in sea level.
One of the first attempts to close the sea level budget compared time series of total sea level from satellite altimetry with the sum of the OBP component from satellite gravimetry and the steric component from Argo floats (Willis et al., 2008). This study showed that between 2003 and 2007 the sum and the total sea level have comparable seasonal and interannual sea level 10 variability, however the 4-year trends did not agree. In that same year Cazenave et al. (2008) found comparable estimates of steric sea level estimated from Argo and from the difference between altimetry and the Gravity Recovery And Climate Experiment (GRACE) observations over 2003-2008. Using the same methods as Willis et al. (2008) the global sea level budget was closed within error by Leuliette and Miller (2009) and Leuliette and Willis (2011) over the period 2004-2008. All of the beforementioned studies used a form of reduced space objective interpolation (Bretherton et al., 1976) to create grids of Argo 15 data. Li et al. (2013) attempted to close the global budget using temperature and salinity grids from Ishii et al. (2006).
While time series of satellite gravimetry and Argo observations became longer and the processing of gravity fields improved, it became possible to look at basin-scale budgets and patterns. Chambers and Willis (2010) compared gravimetry derived maps of Ocean Bottom Pressure (OBP) to those obtained with steric-corrected altimetry, whereas Marcos et al. (2011) investigated the distribution of steric and OBP contributions to sea level changes and looked at basin-scale differences. Purkey (2014) 20 analysed differences between basin-scale OBP from satellite gravimetry and steric corrected altimetry using Conductivity-Temperature-Depth (CTD) profiles over the period 1992-2013. They showed that both methods captured the large-scale OBP change patterns, but that differences occur when deep-steric contributions below 1000 m are not considered. Over the North Atlantic Ocean the OBP trends were found to be statistically equal, but with large error bars for the steric-corrected altimetry results. Von Schuckmann et al. (2014) found global and large-scale regional (a third of the total ocean) consistency in sea 25 level trends of the three systems in the Tropics as long as areas like the Tropical Asian Archipelago are not considered (Von Schuckmann et al., 2014); but they did not manage to close the budget between 30-60 degrees North and argued that the unability of Argo to resolve eddies in the western intensifications caused the difference in trends. GRACE gravity field solutions. Secondly, we address the effect of several processing steps of particularly on gravimetry data in terms of trend, annual amplitude and (residual) time series. For altimetry we briefly discuss the effect of different averaging methods and analyse the effect on the trends of having a latitude dependent intermission bias (Ablain et al., 2015). For GRACE, DDK5-filtered solutions (Kusche, 2007) (Kusche et al., 2009) are compared with the anisotropic Wiener-filtered Klees et al. (2008) solutions. Finally, basin and sub-basin scale budgets are created, problematic areas are identified and potential causes 5 for non-closure are discussed.
This article will describe briefly the data used in Sect. 2. Secondly, the processing of the three datasets is discussed in the methodology section. In Sect. 4 the processed datasets are compared to existing products. The resulting basin and sub-basin scale budgets are described in Sect. 5. In the final section conclusions are drawn based on the results.
2 Data description 10 This section shortly discusses the data from the three observing systems that are used to determine the sea level budgets.
For the determination of the sum of the steric and the mass components of sea level satellite altimetry data are used. The altimetry data are obtained from the Radar Altimetry Database System (RADS) (Scharroo et al., 2012). RADS contains 1 Hz along-track data, which corresponds to an along-track separation of sea level measurements of approximately 6 km. The files contain ranges, orbits and geophysical corrections for all altimeters that have been flown. In this study only the data of the  solutions (Klinger et al., 2016) up to degree and order 90. All three products are provided with full variance-covariance or normal matrices, which allows for statistically optimal filtering. In case of a proper error description, we expect that the results 15 of the CSR 60-and 96-degree solutions give comparable results, except in areas with large gradients in gravity. However, since the variance-covariance matrices are comparable over the periods July 2003-December 2010 and February 2011-July 2013, but the orbit geometry substantially varies, the provided variance-covariance matrices are not expected to give a proper representation of the error, which leads to reduced quality filtering. As a consequence, the months July-October 2004 are excluded from the analysis, when GRACE entered a near 4-day repeat-orbit. The addition of the ITSG solutions, enables us 20 to compare an independent solution computed with a different approach to the standard CSR products. The non-dimensional gravity field coefficients are converted to units of Equivalent Water Height (EWH) before filtering, to make them compatible with the other two observing systems. For comparison, we also used the publicly available DDK-filtered solutions of CSR, however no variance-covariance matrices for those solutions are publicly available. In the processing phase, the Atmosperic and Ocean De-aliasing Level-1B (AOD1B) product is incorporated (Dobslaw et al., 2013), which is based on the Ocean Model 25 for Circulation and Tides (OMCT) and the European Centre for Medium-range Weather Forcecast (ECMWF) model. Monthly averages of the OMCT and the ECMWF are restored after processing to the time-varying gravity field in the form of spherical harmonics (Chambers and Willis, 2010), details on this are described in Sect. 3.3.

Methodology
The data described in the previous section are processed such that they are suited for establishing monthly regional sea level 30 budgets. It implies that the equation  is satisfied within uncertainties, whereh sla is the Mean Sea Level (MSL) anomaly derived from the Jason satellites,h ssla the mean steric sea level anomaly derived from Argo andh obp the mean OBP in terms of EWH derived from GRACE. This section describes therefore the processing strategies for the three observation types from individual measurements to an average over a specified region in the ocean including the propagation of the formal errors.
As far as altimetry is concerned, after computing individual along-track sea level anomalies, two important processing steps 5 are described in this section: a suitable averaging method to come to a time series of MSL for a given area and a way to deal with geographical dependencies of the intermission bias between the two Jason missions (Ablain et al., 2015).
To compute steric sea levels from Argo temperature and salinity measurements the TEOS-10 software is used (Pawlowicz et al., 2012). Since the Argo measurements are non-homogeneously distributed over the ocean, the steric sea levels are first interpolated using an objective mapping procedure to a grid of 1 × 1 degree, before being averaged. 10 Monthly GRACE solutions of CSR and ITSG are provided with calibrated variance-covariance matrices, which allows the use of an anisotropic Wiener filter (Klees et al., 2008). Compared to other existing filters, it strongly reduces the stripes that are still present in the DDK-filtered solutions (as will be shown in Sect. 4, while not reducing the spatial resolution by applying a large Gaussian filter. A fan filter is applied after the optimal filter to reduce ringing artefacts that occur close to Greenland due to the limited number of spherical harmonic coefficients (d/o 60-96).

Jason sea level
Individual sea level anomalies h sla measured with the Jason-1 and Jason-2 satellites are computed with respect to the mean sea surface (mss) DTU13 as:  where a is the satellite altitude, R the Ku-band range and ∆R corr the applied geophysical corrections. The satellite altitude is taken from the GDR-D orbits. The latest versions of the geophysical correction are applied, as listed in Table 2. Sea level anomalies larger than 1 meter are removed from further processing, as in the NOAA GMSL time series (Masters et al., 2012).
In GMSL time series an intermission bias correction is applied, which is determined from the average GMSL difference between Jason-1 and Jason-2 during their tandem phase, in which the satellites orbit the same plane only a minute apart 5 (Nerem et al., 2010). However, the differences reveal a geographical dependence as shown in Fig. 3. Regional sea level budgets established in this study are proner to these geographical differences than when estimating global sea level budgets. This problem is partly corrected for by estimating a polynomial through the intermission differences, which only depends on latitude (Ablain et al., 2015) and is given by: where λ is the latitude and ∆h sla,ib (λ) is the intermission correction. The sea level anomaly h sla,c corrected for intermission differences is then computed as: This correction is only applied to Jason-2 sea level anomalies. The parameters c n , with n = 0, 1.., 4, depend on the applied geosphysical corrections. For the corrections given in Table 2 the values for the parameters are given in Table 1. At halfway  (Nerem, 1995). From here on the latter is referred to as the 'Nerem method'. The gridding method is problematic when using the Jason satellites, because of their large track spacing at the Equator, causing the number of invidual observations per grid cell to decrease at low latitudes (Henry et al., 2014). A solution is to increase the grid cell size, but this has a disadvantage if sea level budgets are constructed over an irregular and/or a small polygon. The latitude weighting method has the disadvantage that it underweights measurements at high latitudes (> 50 degree) (Scharroo, 2006), because it assumes the number of measurements to go to 5 infinity at the inclination of the satellites.
Therefore it is suggested to average the sea level anomalies based on the number of available measurements within a latitude band. The method connects the weights assigned to the measurements to the number of measurements N l in a latitude band l of one degree and the area of the sea surface A l in the following way: 10 These weights are normalized: A MSL anomalyh sla for an area is computed with: whereŵ is the vector of normalized weights andĥ sla,c is the vector of sea level anomalies corrected for intermission differ- 15 ences.
For the error estimation variance-covariance matrices are computed as described in Le Traon et al. (1998). This method separates the long-wavelength errors from the representativity errors due to ocean variability. White measurement noise is not considered here, because it becomes very small when averaged over large areas. Among the long-wavelength errors, we consider the orbit, ocean tide and inverse barometer errors. These errors are assumed to be fully correlated between measurements where λ avg is the average latitude of two measurements. Ultimately, this results in equations for the covariance of respectively measurements in different tracks and on the same track: where ρ ij is correlation computed with the decorrelation time and distance provided above, σ 2 ov is the ocean variability variance and σ 2 lw is the long-wavelength variance. The values for σ ov and σ lw are assumed 100 mm and 15 mm, where the first number comes from typical mesoscale variability (Chelton et al., 2007). By putting these equations in the variance-covariance matrix C sla , the standard errorσ sla for the mean sea level anomaly is computed using: 10 Both the satellite altimetry mean sea level anomalies as well as the EWHs are affected by Glacial Isostatic Adjustment (GIA). In case of altimetric measurements this is corrected by subtracting the geoid trend from an Earth model with VM5a viscosity profile and ICE-6G deglaciation history (Peltier et al., 2015). Errors in the GIA trends are typically assumed to be in the order of 30% of the signal (Von Schuckmann et al., 2014).
Because the CSR gravity fields are created on a monthly basis and the altimetry measurements are averaged over a cycle 15 of approximately ten days, the altimetry measurements are low-pass filtered. A low-pass filter f lp is computed by taking an Inverse Discrete Time Fourier Transform, which results in: with f c the cut-off frequency, which is taken as 12 cyc/year, t is time vector of the altimetry time series and t m the time at the middle of a month. This filter is infinitely long, so therefore we cut it at two months. To obtain a better frequency response the 20 filter is windowed using a Hamming window w H , so that: where L is the length of the window which is two months. The applied filter h lp is then written in the time domain as: 3.2 Argo steric sea level where P 0 is the surface pressure, P is the reference pressure, which is set to 1000 dBar (approximately 1000 meters depth) and g 0 is a constant gravitational acceleration of 9.763 m/s 2 .

5
In the analysis only profiles that reach at least 1000 meters depth are included and at least have a measurement above 30 meters depth, which is the typical depth of the mixed layer. A 'virtual measurement' is created at 1 meter depth, assuming the same salinity and potential temperature values as the highest real measurement, so that the top steric signal is not missed. Only measurements that have error flag '1' (good) or '2' (probably good) are used and the measurements are cleaned by moving a 5x5 degree block to remove measurements more than 3 RMS from the mean. 10 To be able to average measurements monthly over a basin or a polygon, a grid is constructed by statistical interpolation of the steric sea levels at the profile locations based on the method described in Bretherton et al. (1976) and Gaillard et al. (2009). First, a background field is constructed by estimating a model through the 1000 closest measurements of a profile or grid cell location. This model contains a constant, a second-order 2-D lon-lat polynomial and six intra-annual to annual cycles (Roemmich and Gilson, 2009). The background sea surface height is taken as the model evaluated at the grid cell (or profile) 15 location.
Consecutively, the background field vectorĥ ssl,b is subtracted from the sea level estimates, which results in: where δŝ sl are the residuals. The ocean variance σ 2 t is assumed to be 100 cm 2 (typical mesoscale variability (Chelton et al., 2007)), which is close to the average squared RMS-of-fit of the differences of the measurements with the model. These 20 variances are subdivided into three components to represent different correlation scales as follows (Roemmich and Gilson, 2009): which are then used to construct covariance matrices C(d) based on those used for the Scripps fields (Roemmich and Gilson, 2009), such that: and the measurement and representativity error matrix R: where d is a measure for the distance between the profiles p and the grid points g, such that: The parameter a is 1 above 20 degrees latitude and below that it decays linearly to 0.25 at the Equator, in order to represent the zonal elongation of the correlation scale here (Roemmich and Gilson, 2009).
Using the correlations C pg (between profiles and grid points) and C p (between profiles), the weight matrix K is computed 5 as: The weight matrix is then used to compute a vector of steric sea levelsŝ ssl,g for every grid point with in the area: for which also the variance-covariance matrix C ssla , so that: where C g are the correlations between grid locations.
In order to average the steric sea level anomalies, the values are weighted by the cosine of the latitude, which results in: Like for altimetry, the weights are normalized: Eventually these are used to compute the mean steric sea levelh ssl,g and its associated errorσ ssl,g , with and

GRACE ocean bottom pressure
We use the full variance-covariance and normal matrices to filter the spherical harmonic coefficients with an Anisotropic Non-Symmetric (ANS) filter (Klees et al., 2008). This Wiener filter exploits the ratio between the variance of the error and of the signal to filter the coefficients. With the variance-covariance matrices C x and D x , for respectively the errors and the signals, the spherical harmonic coefficientsx are filtered as: 10 Ocean Sci. Discuss., doi: 10.5194/os-2016-50, 2016 Manuscript under review for journal Ocean Sci. Published: 6 July 2016 c Author(s) 2016. CC-BY 3.0 License.
For the filtered coefficientsx of a corresponding variance-covariance matrix C x,of is computed. This is a joint inversion of a static background field, which is set to zero, and the time-varying coefficients, resulting in: The derivation is elaborated in A.

5
The filtered grids contain ringing effects around strong signals over Greenland and the Amazon region, which can have substantial effects on the estimated trends in the ocean, which is discussed in Sect. 4. If averaged over large areas this will have hardly an affect, but on regional scales the ringing should be reduced. To obtain smoother fields, we use a fan filter (Zhang et al., 2009;Siemes et al., 2012), which is given as: 10 where the degree l and order m are related to the maximum degree, l max . For a maximum degree of 60 and 96, this is comparable to a Gaussian filter of 280 km and 110 km, respectively (Siemes et al., 2012). Suppose F f = diag(sinc( l lmax )sinc( m lmax )), then the resulting covariance matrix C x,f f is written as: Note that there is a fundamental difference between filtering the CSR and ITSG solutions. The CSR solutions are computed 15 with respect to a static gravity field, while the ITSG solutions are computed with the respect to a static gravity field including a secular trend and an annual cycle. As a consequence, the CSR spherical harmonic coefficients and signal variance-covariance matrices include the annual and the secular trend while Wiener and fan filtering. The ITSG gravity fields are Wiener-filtered first, then the annual cycle and the secular trend are added back and eventually the fan filter is applied.
Since the degree-1 coefficients are not measured by GRACE, we add those of Swenson et al. (2008) to the CSR solutions.

20
For the ITSG solutions a dedicated degree-1 solution is computed, using the same approach. Furthermore, we replace the C 20 coefficient with satellite laser ranging estimates (Cheng et al., 2013).
The intersatellite accelerations of GRACE are dealiased for high frequency ocean and atmosphere dynamics with the Atmospheric and Ocean De-aliasing Level-1B (AOD1B) product. Monthly averages of the AOD1B are provided as the GAD product for both CSR and ITSG, where the mass changes over land are set to zero. To be able to combine GRACE OBP with inverse 25 barometer corrected altimetry, the GAD products containing the modelled oceanic and atmosphere mass are added back in the form of spherical harmonics. Because the ocean model in the AOD1B product is made mass conserving by adding/removing a thin uniform layer of water to or from the ocean, the degree zero is removed before subtraction from the GAD product to compensate for the mean atmospheric mass change over the ocean, which is not measured by altimetry (Chambers and Willis, 2010).

30
To compute the OBP at a specific grid cell, the 4π-normalized associated Legendre functionsŷ are evaluated at its latitudelongitude location. The OBP of the grid cell h obp is consecutively calculated with: For multiple grid cells the vectorŷ becomes a matrix Y , such thatĥ obp becomes a vector of EWHs. This is written in form: It is possible to compute the grid's variance-covariance matrix C obp as (Swenson and Wahr, 2002): The averaging over an area is equal to that of the Argo grids. Suppose thatŵ are the normalized latitude weights for the grid 5 OBPs, then is the mean ocean bottom pressure in and is its error. 10 Ultimately, the mean GIA OBP over a polygon (degree 2 and higher, as measured by GRACE) are subtracted from the results, to make them compatible with altimetry.

Comparison with existing products
In this section a comparison is made between existing products and the sea levels from altimetry, gravimetry and Argo floats.
First, we compare the MSL time series over the North Atlantic with the existing time series provided by the NOAA Labora-15 tory for Satellite Altimetry (Leuliette and Scharroo, 2010) and we show the effect of a latitude dependent intermission bias.
Secondly, amplitude and trend grids of steric sea level are compared to those computed from Scripps salinity and temperature grids (Roemmich and Gilson, 2009) and the Glorys reanalyses grids (Ferry et al., 2010). Thirdly, the optimally and fan filtered gravimetry grids are compared to the DDK5-filtered gravity fields (Kusche, 2007) (Kusche et al., 2009).

Steric sea level
In Fig. 6 grids for the amplitudes and trends of the steric signal are shown. The Scripps grids (Roemmich and Gilson, 2009) and our solution are solely based on Argo data, while the Mercator reanalyses product Glorys 2V3 assimilates various types of data including altimetry (Ferry et al., 2010), sea surface temperature and Argo. Besides the different input data, the Glorys 10 relies on a volume conserving ocean model, while the other two methods are statistical, the results can differ quite substantially.
Since we use the same correlation structures as Scripps, the resulting grids should resemble each other quite closely. However, to be able to create a variance-covariance matrix between grid cells, it was required to do a 2D-interpolation of the steric sea levels instead of a 3D-interpolation of temperature and salinity profiles. In order to do this, the criteria for removing profiles is, as described in Sect. 3.2 is different for both methods. As a consequence of the 2D-interpolation and the differences in the

Ocean bottom pressure
In Fig. 7 the trends and amplitudes of the CSR DDK5-filtered gravity fields compared with those obtained from the Wienerfiltered ones from CSR and ITSG. Note that the Wiener-filtered solutions are also fan-filtered, as discussed in Sect. 3.3, but will be referred to as Wiener-filtered from here on. Both in the annual amplitude and the trend grids some residual striping To determine how this effects sub-basin scale OBP time series, it is first required to determine the minimum area over which the measurements have to be integrated. GRACE gravity fields have a resolution of typically 250-300 km half-wavelength (Siemes et al., 2012). For small ocean signals after applying filtering procedures, we expect the resolution to be closer to 400-500 km. Argo has approximately one to two floats per 3 × 3 degree box, so its resolution is in the same range as that of GRACE. Jason-1&2 have an inter-track spacing of 315 km at the Equator, which decreases substantially at 60 degrees The month-to-month noise of the CSR Wiener-filtered time series is comparable for all three polygons. The DDK5-filtered time series becomes much noisier for the meridionally oriented polygon, where month-to-month jumps of 10-20 mm occur.
In addition, the DDK5 time series exhibit a substantially different trend in the meridional polygon than the other time series, because the orientation of the polygon is aligned with the residual stripes (Fig. 7). So in terms of trend and noise the DDK5 20 time series strongly depends on the orientation of the polygon. Eventhough the ITSG trend and amplitude grids suffer from striping, they do not become significantly noisier for the meridionally oriented polygon.

Results and discussion
The first objective of this section is to reveal patterns of sea level amplitudes and trends in the North Atlantic and how these resemble for the two different methods: altimetry and Argo+GRACE. Secondly, this section discusses the closure of sea level 25 budgets over polygons of approximately one-tenth of the North Atlantic in terms of trend, annual amplitude and residual variability. It is shown for which regions the budget is closed and possible causes for non-closure are discussed. Thirdly, we focus on the best choice of GRACE filter solutions for the OBP component.

North Atlantic sea level patterns
In Fig. 9 grids of trends and amplitudes computed from Argo+GRACE are overlaid with Jason derived trends and amplitudes 30 at the ground-tracks. In areas where the ground-tracks of altimetry are barely visible, there is a good resemblance between Argo+GRACE and altimetry. The trends from altimetry in the right column of Fig. 9 show a distinct pattern, where positive trends are found below 35 degree latitude and negative trends above this line, with the exception of the North-American coastline. Large trends along the North-American coast are also found by tide gauge studies (Sallenger at al., 2012), where they attribute this to a weakening 5 Atlantic Meridional Overturning Circulation (AMOC). The Wiener-filtered Argo+GRACE solution resemble the trend patterns derived from altimetry measurement better, while the residual stripes in the DDK5 solution are clearly visible. Note that a significantly larger altimetric trend is visible in polygon H, just in front of the Mediterranean. Possible causes will be discussed later.

Sub-basin scale budgets 10
The North Atlantic is split into ten regions, divided in the middle by the Mid-Atlantic ridge, while in the latitude direction trying not to cut through the major oceanographic features, like the salt water tong in front of the Mediterranean and the Gulf

19
Ocean Sci. Discuss., doi: 10.5194/os-2016-50, 2016 Manuscript under review for journal Ocean Sci. Published: 6 July 2016 c Author(s) 2016. CC-BY 3.0 License. Stream, as shown in Fig. 9. As in Sect. 4.3, the size of basins is chosen based on the criterion that the error on the does not exceed 1 mm yr −1 . First, we will discuss three representative examples of time series. Then budget closure in terms of trend and annual amplitude is addressed and the corresponding best gravity filter is determined. Ultimately, the trends, semi-annual and annual signals are reduced from the time series and the best filter choice in terms of residual variability is determined.

Timeseries
Budgets for three representative regions, using the Wiener-filtered OBP solutions, are shown in Fig. 10. The time series for the rest of the regions can be found in the supplementary material. The left plots confirm that the main driver for annual fluctuations in sea level is the steric sea level, but that the trend is strongly influenced by a mass component. On the right side we see that the sum of the components and the total sea level agree to within the error bars, but some problems arise in the Gulf Stream area

Trends
Trends computed from the time series of Fig. 10 are given in Table 2. Close to the Equator (A, B, I and J) and over the whole North Atlantic the trend budget is closed within two standard deviations no matter which of the OBP solutions is used. This confirms the results of Fig. 9. In this region the GIA correction is relatively small and no sharp gradients or strong features are present in the trend grids, which contribute to proper budget closure. Budget closure is also achieved by all filters in the 5 northeast of the Atlantic (F, G). Again, this is a relatively quiet region, with no significant gradients in trends and a small GIA correction. The results of the Argo+GRACE however show are larger spread for the different OBP solutions. It is important to note that especially region F suffers from some ringing artefacts before the fan filter is applied and that the far northeast is not very well covered by Argo floats. The trends of the DDK5-filtered solution in region G are a bit further off than the other solutions, probably resulting from the striping as visible in Fig. 7. 10 In the northwest of the Atlantic the choice of gravity field filter either substantially influences the estimated trends (D and E), or they are just outside of two standard deviations (C) for one ore more solutions. Using the 96-degree CSR Wiener filtered solution, the budget is closed within two standard deviations for all three polygons, but the others do not. For region C, the results of the all filters are resembling quite well, but some are just outside of two standard deviations from altimetry. For region D, the 60-degree CSR Wiener-filtered results are far off, but the other results are close again. In this region, sharp gradients 15 occur not only in the OBP component with the presence of a neighbouring continental shelf, but also in the steric component.
This might lead to leakage of the continental shelf OBP signal or problematic interpolation of the Argo steric sea levels. In addition for both of the beforementioned regions, the GIA correction on the OBP is relatively large. Adding a GIA correction error of 10-20 %, which is smaller than discussed in Sect. 3.3, to the OBP trends would close the budget in these regions for all the solutions, except for the CSR 60-degree solution in region D. In region E, a clear split is visible between the Wiener-filtered 20 CSR solutions, which close the budget and the other two solutions, which do not close the budget. The difference in results could be caused by the filter not being able to handle the large gradients (Klees et al., 2008) in OBP within this region (Fig. 9).
However, if we would again add only a 10-20 % GIA correction error, it would suffice to close the budget for all filters.
Ultimately, only the budget in region H cannot be closed with any of the solutions and there is no strong GIA signal present, which could be responsible for a large bias. In addition, the sea level in this polygon does not exhibit any strong gradients and 25 the number of Argo floats is substantial. This excludes interpolation or filtering problems. Therefore, we argue that this can be explained by a deep-steric effect, that could be related to variations in the export of saline water from the Mediterranean (Ivanovic et al., 2014), which is not captured by Argo.
In conclusion, it is possible to close the sea level budget within two standard deviations for nine-out-of-ten regions using the Wiener-filtered 96-degree solutions. If a 10-20 % GIA correction error is taken into account, the budget for nine-out-of-30 ten polygons is also closed for the DDK5-filtered solution and the Wiener-filtered ITSG solution. This also suggests that the commonly assumed GIA correction error of 20-30 % (Von Schuckmann et al., 2014) is probably overestimated. Table 2. Trends of total sea level (mm yr −1 ) and their standard deviations from altimetry (Jason) and the sum of steric and OBP from Argo (A.) and GRACE (CSR, ITSG) for different filter solutions. NA is the trend for the complete North Atlantic between 0-65 degree latitude. A 0.4 mm y − 1 drift error is taken into account for altimetry based on the comparisons with tide gauges (Mitchum, 1998(Mitchum, , 2000.

Annual signal
We indicated that the seasonal cycles are primarily caused by steric variations in sea level (Fig. 10). By comparing the first column with the last column in Table 3, it becomes clear that in most cases an additional OBP signal is required to close the budget in terms of annual amplitude. The discrepancy between Argo and altimetry for the whole North Atlantic reveals that on average in-phase OBP signals with an amplitude of approximately 7 mm are required to close the budgets, which is in line 5 with the modelled results of Tamisiea et al. (2010). They modelled, using fingerprints, amplitudes of OBP ranging from 3-12 mm, and phases (not shown here) between day 210-330, which is in-phase with the steric signal. amplitude with respect to altimetry in region D. In region B the estimate of ITSG+Argo is relatively small and in region D the CSR DDK5+Argo also relatively large. Note that the number of Argo floats in region B is often small (Fig. 1) and that

22
Ocean Sci. Discuss., doi:10.5194/os-2016Discuss., doi:10.5194/os- -50, 2016 Manuscript under review for journal Ocean Sci. Published: 6 July 2016 c Author(s) 2016. CC-BY 3.0 License. large gradients in the steric sea level in region D could cause interpolation problems for steric sea level. Secondly, in both northern polygons E and F both combinations of Argo+GRACE underestimate the amplitude compared to altimetry. Why this underestimation occurs is not completely clear. A likely culprit is the gravity field filtering, but yearly deep convection events in these regions (Våge et al. , 2009), which transport surface water to depth below 1000 m, and the limited number of Argo floats, could also be contributing factors.

5
Using the ITSG solutions it is also possible to close the budget on the scale of the whole North Atlantic (last row of Table 3).
The Argo+ITSG performs best in terms of amplitude budget closure in most regions, eventhough often characterized by sligthly smaller amplitudes than those derived from altimetry. This suggests that there is either a long-wavelength underestimation of the amplitude in GRACE, an over estimation in altimetry, or a missing steric effect in Argo. This is in line with Storto et al. (2015), where on a global scale, steric sea levels computed from reanalyses and gridded T/S fields are found to be smaller than

Residual variability
Time series for the same regions as in Fig. 10 are shown in Fig. 11, but their trend, semi-annual and annual signals have been This suggest a link between the latitude of the Gulf Stream and sea level temperatures in the East of the Atlantic. In region B, 10 we also observe a small interannual effect by altimetry and Argo. However, the amplitude of the signal is larger for altimetry than is captured by Argo, which suggests either some interpolation issues in an area without many Argo floats or a deep-steric effect.
Using any of the filtered CSR or ITSG solutions, it is possible to detect the interannual variability described, probably because most of the signal is of steric origin. However, for the interannual signals that are less pronounced, or for high frequency  behaviour of sea level there are some differences between the OBP solutions. Table 4 shows the fraction of variance of the residual signal of altimetry (trend, semi-annual and annual cycles removed) explained by Argo+GRACE.
The third column indicates that Argo in combination with Wiener-filtered 96-degree CSR solution does not explain much of the residual variance, but mostly introduces additional noise, which causes the negative values. Using the DDK5-filtered OBP the explained variance increases, but the best performance is obtained with the Wiener-filted 60-degree CSR and especially the 5 90-degree ITSG gravity fields. The last column shows that after reducing the trend, and the semi-annual and annual signals, between 24-53 % of the residual signal can be explained by the combination of Argo and ITSG. It is remarkable that for the whole North Atlantic (last row), no variance is explained by the Argo+GRACE, primarily due to the absence of a clear interannual signal. closure in nine-out-of-ten regions. In the considered regions also the DDK5-filtered CSR solutions and the Wiener-filtered ITSG solution appear to close just as many budgets when a 10-20 % GIA correction error is added. The results of the DDK5 filter however, strongly depend on the orientation of averaging area due to residual meridional striping. The strong resemblence between trends also indicates that the errors on the GIA model are probably smaller than the commonly assumed 20-30 %.
Furthermore, a large difference in trend between altimetry and Argo+GRACE is observed in front of the Mediterranean where 5 only a small GIA correction is applied. We believe that this originates from steric effects below the considered 1000 m, where saline water enters the Atlantic from the Strait of Gibraltar and dives to large depths. Further research is needed to confirm this hypothesis.
The 60-and 96-degree Wiener-filtered CSR solutions appear to underestimate the amplitude of the annual signal substantially. They also suffers from what appears to be leakage around the Amazon and Sahel, regions with a substantial annual 10 hydrological cycle. Using the DDK5-filtered gravity fields and the ITSG solutions, the sum of the steric and OBP components becomes significantly closer to that of altimetry, with closure in seven-out-of-ten regions. However, it must be noted that the altimetry signals tend to be slightly larger. This is likely due to partly destruction of the signal by filtering of the gravity fields or limited Argo coverage, or in some regions deep-steric signals.
By removing the semi-annual and annual signals and trends interannual variability can be detected. Since most of the 15 interannual variability in the North Atlantic is contained in the steric component, the type of filter on the gravity fields is not really important. However, if we look at differences on a month-to-month basis, high-frequency variations or small interannual fluctuations in OBP, it is best to use the 60-degree anisotropic Wiener-filtered solutions or the ITSG solutions, because the fraction of explained variance of the altimetric sea level time series by the sum the components using these solutions is largest.
Using the ITSG solution, 24-53 % of the variability in the altimetry-derived sea level time series is explained. The 96-degree 20 Wiener-filtered CSR solution only introduces noise and explains virtually no residual variability of the altimetry time series.
Especially in the 4-day repeat-orbits in 2004 and even the months around them, the Wiener-filtered solutions do not give proper estimates of the OBP components, which partly contributes to a lower explained variance.
To summarize, using the ITSG Wiener-filtered solution the trend budgets close when a GIA correction of 10-20 % is applied.
They perform, together with the standard DDK5-filtered CSR solution, best in terms of annual amplitude budget closure. 25 Additionally, the combination of ITSG OBP and Argo steric sea levels explains the largest fraction of variance in altimetry time series. Based on this, the best option to establish budgets, at scales considered in this paper, is the Wiener-filtered ITSG solution. However, due to residual striping in the trend grids from the 'static' background field that are added back after Wienerfiltering, one must take care when averaging OBP over even smaller regions, or meridionally oriented polygons, which is a even a bigger problem for the standard DDK5-filtered CSR solutions.

Appendix A
The Wiener filter is in principle a joint inversion between the spherical harmonic coefficients of the background fieldx b and those of the time-varying gravity fieldx. Suppose that C x is the error variance-covariance matrix ofx and D x the signal variance-covariance matrix, then the filtered coefficientsx f are expressed as: Assuming the spherical harmonic coefficients of the background field are zero, this equation reduces to Eq. 27. Its filtered variance-covariance matrix C x,f is computed using: Since the matrices (C −1 x + D −1 x ) −1 and C −1 x are symmetric, it is possible to simply change the order underneath the transpose sign and leave: which is further simplified by using the identity C −1 x C x = I to: 10 Finally, this equation is rewritten, such that: which eventually is simplified to Eq. 28.
the service of the employees of Mercator for delivering the Glorys reanalysis product. We wish to express our gratitude to the RADS teams for constantly maintaining and updating their database. This study is funded by the Netherlands Organisation for Scientific Research (NWO) through VIDI grant 864.12.012 (Multi-Scale Sea Level (MuSSeL)).