Weighing the ocean with bottom-pressure sensors: robustness of the ocean mass annual cycle estimate

We use ocean bottom-pressure measurements from 17 tropical sites to determine the annual cycle of ocean mass. We show that such a calculation is robust, and use three methods to estimate errors in the mass determination. Our final best estimate, using data from the best sites and two ocean models, is that the annual cycle has an amplitude of 0.85 mbar (equivalent to 8.4 mm of sea level, or 3100 Gt of water), with a 95 % chance of lying within the range 0.61– 1.17 mbar. The time of the peak in ocean mass is 10 October, with 95 % chance of occurring between 21 September and 25 October. The simultaneous fitting of annual ocean mass also improves the fitting of bottom-pressure instrument drift.


Introduction
The total mass of water in the oceans fluctuates with seasonal changes in continental water storage.A measure of the annual cycle of water exchange is of widespread interest and has been estimated in at least nine studies using data and models including satellite gravity, hydrology, ocean steric height and satellite altimetry.Amplitude estimates have a wide range of 5.5-9.4mm (Vinogradov et al., 2008;Wouters et al., 2011), and phases from 259-303 • (Siegismund et al., 2011;Rietbroek et al., 2009).(A phase of zero represents an annual peak at the start of the year, so these correspond to a maximum ocean mass between 19 September and 3 November.) It is theoretically possible to use a single bottom-pressure sensor to monitor changes in the total mass of water in the oceans independently of satellite gravity measurements and, if the sensor location is carefully selected, with little dependence on hydrology models (Hughes et al., 2012).There is also a great deal of interest in monitoring interannual or decadal variations and long-term trends in ocean mass, but existing bottom-pressure sensor technology makes this extremely difficult, as the instruments, usually deployed for high-frequency tasks such as tsunami monitoring, suffer from non-linear drift of the order of centimetres per year.The drift can vary even between redeployment of the same instrument (Watts and Kontoyiannis, 1990;Polster et al., 2009).However, a good determination of the annual cycle of ocean mass change is still valuable for constraining models of hydrology and ocean dynamics, and in providing an independent measurement for comparison with GRACE and altimetry.Conversely, knowledge of the annual cycle also allows us to improve estimates of the non-linear contribution to bottom-pressure instrument drift.
Based on a pair of bottom-pressure sensors moored in the Pacific, Hughes et al. (2012) estimated an amplitude 8.5 mm (equivalent to a global average pressure of 0.86 mbar, or about 3200 Gt of water) and phase 262 • (22 September).This lies within the envelope of results from other studies, but no formal attempt was made to put error bounds on this number, although it was noted that a similar value could be derived using different ocean models.
Ocean dynamics aside, the bottom-pressure cycle measured at a given site is affected by the crustal deformation and gravitational effects caused by the mass change on land.The bottom pressure needs adjusting everywhere to derive the global ocean mass change, but at certain sensor locations in the central Pacific the result is uniformly biased high.At these locations the effect on local bottom pressure is almost independent of the origin of the additional water mass, which Published by Copernicus Publications on behalf of the European Geosciences Union.could equally come from Greenland, Antarctica or the Amazon.This was the reason for the choice of sites S and N of Hughes et al. (2012), as the nearest existing data to the optimal region.
In this paper, we will apply the Hughes et al. ( 2012) "weighing" technique to existing bottom-pressure data from 17 moorings (including sites S and N) at tropical sites around the world.Many of these represent sites that would not be considered ideal, either due to increased variance of the dynamic bottom-pressure signal or dependence upon the specifics of the continental hydrology.We use ocean models to remove local dynamics, including self-attraction and loading corrections.We use long-period tide models and atmospheric pressure data to correct for those components omitted from ocean models.
We will first, independently of any hydrology model, show the range of annual cycles in local pressure that arise from mass exchange, and how the error bounds vary between sites according to the quality of available data.Then, with an hydrology-and-atmosphere model based on GLDAS (the Global Land Data Assimilation System, Rodell et al., 2004), we account for the expected differences between sites and make use of data from these sub-optimal locations, reducing error bounds on the predicted annual cycle of global ocean mass.We will describe three techniques to derive error estimates, and how to combine data from multiple sites.Our best estimate is that the ocean annual cycle has an amplitude of 0.85 mbar and phase of 280 • (10 October), with 95 % of results within 0.61-1.17mbar or 261-294 • (21 September-25 October).
2 Calculating annual cycle in bottom pressure at the 17 sites individually: method

Contributions to bottom pressure
The bottom pressure p rec measured by a given sensor can be decomposed as where p drift is sensor drift, p dyn is the change due to ocean dynamics, p t is due to tides, p a is the atmospheric pressure averaged over the ocean and p m is due to the change in ocean mass due to precipitation, evaporation, grounded ice melt and river runoff.Some function F s relates the global-average mass change m o to the pressure felt at an individual site.In general, F s , an adjustment due to the changing geoid and crustal loading, is dependent not just on the site location but on the distribution of ocean and continental water mass, and cannot be assumed to be stationary.However, for the annual component we may assume that such a uniquely invertible function exists for a given site, (2) In this part of the paper we seek only to determine the annual cycle of the local pressure p m at each site, and not m o itself.This annual cycle of p m is calculated completely independently of any hydrological or atmospheric model, or GRACE data.In Sects. 4 and 5 we will employ a hydrological and atmospheric model to determine F s at annual timescales for each site -a rather small correction, at certain sites -and hence the annual cycle in ocean mass.

Sensor locations
For this paper, we use data from the US National Tsunami Hazard Mitigation Program (González et al., 2005), downloaded from http://www.ndbc.noaa.gov/dart.shtml.These were, at the time of selection, all the sites available in the open ocean within 18 • of the equator with records from 2001 onwards.Where records have been converted to equivalent metres of sea level, we revert them to pressure in mbar using the original conversion factor.There are 17 sites, at locations shown in Fig. 1 and Table A1.Sites 17 and 11 correspond to sites N and S of Hughes et al. (2012).These two sites are the earliest, starting in September 2001 and January 2003 respectively; the other sites have deployments between 2006 and 2011.At most sites there are multiple instrument deployments.The deployments varying in length from a few months to over 2 years, and there are many gaps in data both between and during deployments.Sites 3, 8, 9, and 16 have particularly short records.
We initially do an approximate detrending, in order to fit and remove tides with daily or shorter periods and the fortnightly tide components M f and M sf with periods 13.66 and 14.77 days.Then we replace trends ready to calculate a more precise drift fit in combination with the annual signal, as described in Sect.2.2.The data with fortnightly and shorter tides removed is shown in Fig. 2.

Long-period tides
There remain other tidal constituents with monthly or longer periods in p t which can be removed by modelling rather than fitting, and thus avoiding contamination by the sensor drift.(In fact we tried both fitting and modelling for removing fortnightly and monthly equilibrium tides, and found less than 0.01 mbar difference in results for the annual cycle.) The first of the long-period tides we consider is the Pole Tide η p , which arises from variations in the apparent location of the Earth's rotational axis, and hence variations in the centrifugal potential.It has main periods of 12 and 14 months and amplitudes of up to 0.7 mbar at the bottompressure recorder sites.We calculate this using the derivation given by Desai (2002).Over the ocean, .We subtract a constant and linear trend from p x , p y , in order to focus on shorter period signals.
There are also long-period equilibrium tides with amplitude of up to 1.5 mbar over the recorder sites.These are of the form where and the coefficient z 0 corresponds to the components with periods 18.6 years, annual, semi-annual and 4-monthly, and z 1 , z 1a to monthly components.The time series for z 0 , z 1 and z 1a are synthesised using components tabulated by Cartwright and Tayler (1971); Cartwright and Edden (1973).

Self-attraction and loading of long-period tides
The presence of additional water causes crustal deformation and geoid changes that change the ocean depth, both locally and globally.Generically, loading effects such as these are sometimes referred to as self-attraction and loading (SAL).So as well as the effect of the external changes in gravitational potential, we account for SAL corrections to the longperiod tides.
We do this following the method using Green's functions described by Stepanov and Hughes (2004).To calculate E for each η we apply the technique described by Agnew and Farrell (1978), who applied it to an input potential equivalent to our η c .The ocean at a point (φ a , λ a ) feels the effect of the sea-level change at each other point, η(φ b , λ b ), according to a function G(α) where α is the angular distance between a and b.We assume that the additional mass at the point distributed in a circular cosine bell hump, i.e. mass is proportional to (1 + cos(πr/r x ))/2 where r is the arc length of α and r x = 0.25 • is chosen as the grid size.G has three components, G = G 1 + G k − G h , which are the Green functions for vertical seafloor displacement due to loading, G h ; vertical geoid displacement due to indirect attraction, G k ; and vertical geoid displacement due to direct attraction of the load, G 1 .These are taken from interpolations of the Green functions given by Francis and Mazzega (1990).
We integrate the attraction over all points b to give the global function E(η), the additional mass at any point due to the input sea level η: Then, since the correction E(η) itself changes the geoid, the process must be iterated.Eventually we converge on η tot (φ, λ), an equilibrium response to the forcing by the spherical harmonic η, which is self-consistent under loading and self-attraction, so η tot ≈ E(η tot ) = E(η).This is done for η = η x , η y , η c and η 3 , and the result is consistent with the maps of Agnew and Farrell (1978) for η c and Desai (2002) for η x and η y .The process only slightly changes the analytical functions calculated without applying the SAL correction, in most areas scaling them by about 1.25.

Ocean dynamics and atmospheric pressure
To remove the local ocean dynamics p dyn from the bottompressure measurements, we use ocean models.For this study we have available five models, as detailed in Table 1.ECCO and the 1/12 • NEMO runs provide the best overlap with the bottom-pressure data, and results below are based on ECCO unless stated otherwise.
We subtract the global spatial average of bottom pressure at each time from the model data, to remove any added mass, as well as removing artifacts due to the model's Boussinesq approximation.
Figure A1 shows the annual cycles of local bottom pressure that are found in the models at each site, with most sites  (Blaker et al., 2014).
having an amplitude of about 0.5-1 mbar.At some sites (e.g.11) the models give consistent annual cycles, but at site 2 there is almost a factor of 2 difference in amplitude between ECCO and NEMO 1/12 • and at site 17 there is a range of 76 • between OCCAM 1/12 • and NEMO 1/12 • .This is one of the greatest areas of uncertainty in this study.
The atmospheric pressure averaged over the global ocean, p a , needs to be removed from the bottom-pressure record, and this is done using the ECMWF analysis data set provided as a satellite altimeter product by Aviso.p a has an annual amplitude of 0.61 mbar and phase 186 • , peaking at 7 July.

Self-attraction and loading of dynamic pressure
The ocean models assume a constant gravitational field and ocean floor, but in the real world water masses cause crustal deformation and changes to the geoid as described above.For the p dyn component, we correct the model data for these SAL effects, in a similar way to the calculation by Tamisiea et al. (2010).The SAL correction to the dynamic pressure from model data has an annual amplitude of up to 0.2 mbar at these sites.
This correction was not made in the previous study (Hughes et al., 2012) and accounts for a large part of the phase difference between that study and this one.

Removing non-linear sensor drift
Currently available bottom-pressure sensors suffer from drift that can be larger than the annual cycle of bottom pressure (Fig. 2).We follow Watts and Kontoyiannis (1990) and Polster et al. (2009) in assuming that the drift is an initial decaying exponential and a long-term linear, that is of the form where τ is the time in days since the start of the deployment.To prevent our fitted drift absorbing the annual cycle (we cannot prevent it absorbing the long-term trend without further information on instrument drift), we use an iterative process as follows.We write p m = p ANN m + p noise , i.e. the annual signal we seek plus some signal which includes high frequency ocean mass changes and errors in the ocean model.Rearranging Eq. ( 1) as we guess p ANN m (we start with amplitude = 1 mbar, phase = 0) and find coefficients to fit p drift to p res .If we have correctly guessed p ANN m , then there will be no annual signal in the dedrifted residual p noise .So we fit another annual p ANN1 m to p noise , and use this to adjust our estimate of p ANN m .We iterate with fresh attempts at p ANN m and the drift fit, until there is no annual signal left in p noise , and p ANN m has converged.Convergence to a tolerance of 0.005 mbar (amplitude of adjustment between iterations) is usually achieved within about 10 iterations; "no convergence" is declared after 80 iterations.
For most sites in our test set, there is more than one deployment, although there may be gaps between and during deployments.In these cases the iterative procedure involves all deployments at a site, simultaneously fitting individual drifts and a single annual cycle.An example of the fitted drift is shown in Fig. 3.
3 Calculating annual cycle in bottom pressure at the 17 sites individually: results

Requirement for iterative fit for dedrifting
Figure 4 shows the fitted annual cycle of bottom pressure at each site when the recorder drift is calculated (a) with the iterative procedure outlined above and (b) with a least-squares fit of an exponential-plus-linear to the raw recorder data p rec .At some sites (1, 2, 4, 10, 11 and 12) the iterative fitting makes fairly small differences, but the fitted annuals at sites 5, 13, 14 and 15 are changed substantially.
Site 15 provides a particularly clear example (see Fig. 3).The apparent decrease of the bottom-pressure record in late 2008 is due to the coincidence of the annual ocean mass decrease.The drift fitted to the raw data is 21.54 − 0.05τ − 47.47e −τ/84.33 , which is decreasing at the end of the record.The iterative fit allows for the sinusoidal contribution of the ocean mass and other variables, and the resulting drift is −3.76 + 0.05τ − 19.41e −τ/37.24, increasing throughout the record.Using the raw-data drift fit would result in a bottompressure error of over 3 mbar for the end of the record.With a differing sign for the linear part, serious error could result from any extrapolation of the raw-data drift fit.This removal of a real pressure signal occurs for the first deployment of sites 13, 14, and 15, all of which are only 16 months long and finish in April 2008 at the minimum of the ocean mass cycle, maximising the risk of the annual signal contaminating the drift fit to the raw data.For the first deployment at site 13 the raw fit is over 6 mbar too low at the end of the record.
It is worth remembering that we have not attempted to distinguish between the bottom-pressure recorder drift and any long-term trends in ocean mass.So this fit represents a combination of recorder drift and any trend or variability in ocean mass longer than 1 year.Trends in other components of Eq. ( 1) may also remain.

Range of results across test sites
The amplitude and phase of the annual cycles in p ANN m for each of the test sites is shown in Fig. 4a, plotted with a phase of zero at the top representing an annual peak at the start of the year.Site 11 (S) has amplitude 0.93 mbar and phase 272 • peaking at 3 October, slightly larger and later than the result of Hughes et al. (2012) for the combination of sites S and N.For sites 3, 8, 9 and 16, all of which had very short deployments, the iteration did not converge.
At this stage, some of the results differ wildly from other measurements of the ocean annual cycle, with a scatter of 6 months of phase and several times amplitude.They are not expected to be exactly the same as these are measurements of residual pressure cycle (p ANN m ) at each site, with geographic variation in ocean mass still included, i.e. the inverse of the function F s still to be applied to convert to global average ocean mass m o .But it is clear that the scatter is larger than that predicted by the ocean mass model (Fig. 4c).It is perhaps not surprising that the short deployment at sites 4 and 5 produced implausible results, but more so that there is so much difference between neighbouring sites 13 and 15.However, we have yet to examine the error bounds on these measurements and as we will see in following sections, the bounds on some sites are much larger than on others.

Sensitivity to noise
If all dynamical signals were perfectly modelled and removed, then it would always be possible to distinguish annual cycles from sensor drift.However, the presence of noise means that our ability to distinguish these two signals will depend in a complicated way on the type and amplitude of the noise, the nature of the drift, and lengths of time series.
To test this, we produce a random noise signal p noise1 with similar frequency spectrum to the residual p res .To produce simulated time series whose stochastic properties closely match the real data we first form the power spectra from the 17 series.Given that the series contains many gaps we use a combination of the Lomb-Scargle periodogram (Lomb, 1976;Scargle, 1982) and Fourier spectral analysis on segments that were over 5000 epochs (52 days) long to produce a good representation of the spectra.We find that the spectra follow a power law shape, that is the power is proportional to f α , for frequency f and spectral index α.However, the slope appears to change at around a period of 2 days and at the lowest frequencies the spectra is essentially flat, probably due to removing the drift from the series.We estimate the low-frequency spectral index (period greater than 2 days) by averaging the series to daily estimates and using maximum likelihood estimation (MLE).All 17 estimates of the spectral index are then averaged to come out with a spectral index of −1.3.We cannot estimate the high-frequency spectral index using MLE because the sub-daily part of the spectrum is biased by remaining power at tidal frequencies.The spectral index is therefore estimated from the power spectra to be around −1.8.We simulate the noise using the discrete simulation method described by Kasdin (1995) where the impulse response function is created to mimic the −1.8 spectral index at high frequencies and −1.3 index at periods greater than 2 days.Finally, to recreate the flattening at low frequencies we also remove trends from the simulated data for randomly sized segments equivalent to that seen in the real data.
We do not create any extra power at tidal frequencies.
We then adjust each of the 17 original bottom-pressure signals using p rec − p noise + p noise1 -the same p noise1 for every site -and redo the iterative drift and annual fit.This is repeated for 100 noise signals to give a spread of results for each site.Although the models perform better at some sites than others, the same noise spectrum is used for each site.This enables us to directly compare the spread of results at one site with another.
Figure 5 shows the scatter of the resulting annual fits for each site.For sites 3, 8, 9 and 16 only a few results are plotted, this is not necessarily because the amplitude is greater than the range shown but because the iterative fit diverged.The fits with a close grouping, such as 1, 2, and 11, are those for which the annual fit has been robust to the specific noise estimate.Some of the surprising results in Fig. 4a, such as sites 4 and 15, are those for which there is a large spread here.

Comparison of ocean models
The models for ocean dynamics, NEMO, ECCO and OC-CAM, have similar annual cycles of model bottom pressure (p dyn ) for some sites (e.g.sites 7, 10, 11) but agree less well at others (e.g.sites 2, 15) (see Fig. A1).This increases the error bound on the ocean mass for the latter sites, suggesting that more weight should be given to those sites with good model agreement.To calculate this rigorously, we repeat the noise calculation in Sect.3.3 for all models.These have the distributions as shown in Fig. 6.The OCCAM and NEMO4 models have shorter overlaps with the data, and hence provide shorter series for the annual fitting, which may account for the poorer performance of these models.The OCCAM models are only used at sites 11 and 17.NEMO4 has a short overlap with the data at most sites, but overlaps more than 16 months only at sites 1, 11 and 17.At sites 4 and 5 there are only 14 months of data, insufficient for a reliable estimate of the annual cycle.Site 16 has a long gap between deployments, with only 15 months of data in total.We see that the result at site 1 is consistent between models.In Sect. 3 we focussed on p m , the annual that can be retrieved from bottom-pressure records once other known signals had been removed: atmospheric pressure, long-period tides and ocean dynamics.Of these, the last is probably the least known and would introduce the greatest spatial variability into the bottom-pressure records.However, even once these "known" signals are removed, hopefully leaving just the pressure signal due to mass flux, we would not expect p m to be globally uniform (e.g.Clarke et al., 2005).During the annual water cycle, some continental locations store large masses of water, and the SAL effects of this additional mass should be accounted for.Ideally, one might locate bottom-pressure sensors far from land, so that the location of the variation of water mass on the continents would be irrelevant.Hughes et al. (2012) demonstrated that no matter the location of the mass change on the continents, areas of the Pacific Ocean experience nearly the same change in bottom pressure.They divided the continental areas into 2283 separate regions, and looked at the change in bottom pressure associated with a water loss from each region that would be associated with a one millimetre globally averaged sea level rise.While areas near the mass loss experience a bottom-pressure decrease, a section of the Pacific experiences an increase only ranging between 0.9 mm and 1.3 mm whatever the locations of mass loss on the continents.On average, this region experiences a higher-thanaverage sea level change of 1.15-1.18mm.
Unfortunately, the bottom-pressure records are not located in this ideal area, as coastal locations are more relevant in their role in the tsunami warning system.Thus, we expect that crustal deformation and geoid changes can introduce spatial variations into the amplitude and phase of the annual mass signal.We will now introduce how we account for this when trying to find a globally averaged value.

Modelling the spatial variability of p m
The change in bottom pressure p m due to the annual cycle of mass change m o , depends upon the location of the bottom-pressure sensor and pattern of mass change on the continents, according to some function F s .Similar to Hughes et al. (2012), we follow the approach of Tamisiea et al. (2010), using hydrological and atmospheric models to estimate F s as where p h is the model prediction of pressure for that site, and m h is the model estimate of m o .The dominant contribution to the ocean mass comes from the change in water storage on the continents.However, the atmosphere also contributes to the global annual cycle in ocean mass, both by storing water mass and by pressure changes introducing loading changes on the continent.In this calculation, we assume that the Earth responds elastically, which is reasonable for the annual cycle.
For our primary estimate of this effect, we use the GLDAS/Noah data version 1 (GLDAS-1) (Rodell et al., 2004) over the period December 2001 to November 2010.Since this does not extend to Antarctica, and the data in Greenland should not be used, we add the same component from GRACE for these regions as used in Tamisiea et al. (2010) andVelicogna (2009), although this has a rather small effect.For the atmospheric data, we use the monthly-mean surface pressure fields from the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalyses (Kalnay et al., 1996) over the same period.The resulting annual cycle estimate at each of the 17 bottom-pressure locations, as well as the global average, is shown in Fig. 4c.Compared to the global average, 0.86 mbar and a phase of 265 • (25 September), most of the sites show a larger amplitude.Many of the sites (4,(11)(12)(13)(14)(16)(17) show a similar amplification of amplitude with rather small phase change.However, other locations indicate quite different behaviour.Most notable is site 2, located in the Indian Ocean.The large water storage on the continents in the region, with a not-too-different phase with respect to the ocean's maximum, leads to a much larger amplitude (1.21 mbar) and later phase (277 • , 8 October) at this site.The actual values displayed in Fig. 4c are not as important as the scaling relationship one can infer between the globally averaged annual mass variation and the changes at each bottom-pressure recorder location.As long as the distribution of water is correct, both spatially and temporally, then the scaling inferred from the results should be independent of the actual value of the global average.This is the critical assumption employed here.We assume that the relationship between local pressure change and the global average is the same for the model and the observations.

Correcting for spatial variability of p m
The values of p h for each site and m h provided by the hydrology and atmosphere model are shown in Fig. 4c.Then our estimates of m o from each site are given by The calculation is done using complex variables to treat the amplitude and phase of the annual signal together.This In order to quantify the noise distributions using kernel density estimators on the sine and cosine coefficients of the annual signal, we use the function kde2d submitted to the Matlab file exchange by Zdravko Botev (Botev et al., 2010).The distributions for sites 4 and 11 in Fig. 7 illustrate how the results for some sites are much more susceptible to noise than others.Sites with very few converged results are omitted, but note that the probability density functions (PDFs) for these sites will be close to uniform, so will not substantially affect later results.After correction for the spatial variability of p m , sites 2, 6, 7, and 11-14 all show a focussed peak in results with an annual of around 0.9 mbar and phase in early October.Only sites 1 and 10 show a focussed peak with results elsewhere.Site 10 has only one good deployment and the results are slightly smeared.The high amplitude of site 1 is not easy to explain, but notice that its location in the Caribbean Sea is not ideal.We speculate that it may be subject to local ocean dynamics poorly captured by these models, although ECCO and NEMO12 are in reasonable agreement.
The noise spread of the results is much greater than the correction for spatial variability, and it seems unlikely from these results that we could reverse the calculation to detect meaningful spatial variation in continental water storage from the bottom-pressure record.

Combining all the sites
We can combine the sites to make use of as much information as possible, by first using the hydrology and atmosphere model to adjust from p ANN m to m o .Since the PDFs for each site are independent realisations of ocean mass cycle estimates, the PDF for a combined estimate is given by the product of these PDFs, renormalised to integrate to 1.In effect, this gives more weight to those sites with a narrower noise spread, in which we have greater confidence.The resulting PDF is shown in Fig. 8a.The peak has an annual amplitude of 0.93 mbar (3400 Gt) and phase of 287 • (18 October), > 95 % of results fall inside an amplitude range of 0.79-1.06mbar and > 95 % within a phase range of 266-293 • (26 September-23 October).
Including the NEMO12 model as well as ECCO leads to a result that is smaller and earlier than the result for ECCO alone, with peak amplitude 0.82 mbar and phase 282 • (13 October).The smaller amplitude appears to be due to the difference between NEMO12 and ECCO at site 2 (see Fig. 6).

Comparison of different hydrology models
Thus far we have described results using the GLDAS-1 data for the hydrology model to derive the function F s .We have also tested the GLDAS-2.0 data (also plotted on Fig. 4c), which uses an updated version of the NOAH model and, more importantly, different meteorological forcing.While the whole GLDAS-1 model time series is forced by a mix of meteorological data sets, over the period of this study the forcing is consistent and includes high-quality observational precipitation and solar radiation.GLDAS-2 uses the Princeton meteorological forcing data set, which is a bias-corrected reanalysis product.
In this case the global-average ocean mass predicted by the combined hydrology and atmospheric data sets, m h , is very different, with an amplitude of only 0.52 mbar (peaking 24 September).(The ECCO ocean model is used here.)But as seen in Fig. 9 the relationship between m h and the local pressure p h is similar, so F s is little changed at most sites.The largest difference is at site 2 in the Bay of Bengal, for which GLDAS-2.0 predicts a much larger amplitude F s than GLDAS-1.When the results for different sites are combined, the most probable prediction for the annual of ocean mass is 0.84 mbar with phase 277 • , 7 October, earlier than the result with GLDAS-1 and with a smaller amplitude.The spread of results is similar.We believe this is largely because of the increased shift of site 2, which brings it closer to the centre of the spread of other sites' results than with GLDAS-1.
5 Calculating ocean mass from all deployments at the 17 sites simultaneously

Simultaneous fit
As there is only one ocean mass cycle to be determined, it makes sense to consider the alternative approach of calculating a simultaneous fit to all deployments, rather than treating each site separately.This also allows us to improve the fit to  drifts at sites for which all deployments are short, as the annual cycle used in the fitting process will be constrained by longer records from other sites.
To combine records we apply the procedure outlined in Sect. 2 with a single guess at m o , the annual cycle in ocean mass.We use F s (from GLDAS-1) to give p ANN m , the annual cycle in bottom pressure at each site.(Observe that this requires us to have F s , the relative scaling on each location at this stage, although not an absolute value for the ocean mass.)We fit linear-plus-exponential drifts to p res for each deployment, to give p noise , then use a least-squares fit on p noise for all deployments to find any annual in the residuals.This is used to adjust our guess at m o and we iterate to minimise the annual signal remaining in p noise .
The multiple-site drift fitting makes only a small difference to many deployments, but it does improve the fit to sites with very short deployments (e.g. 3, 9).On the first deployment (2007)(2008) of site 15 it has the effect of reducing the effect of the annual signal and shortening the timescale of the exponential part (see Fig. 3).The drift there changes from 1238−0.50τ−1251e −τ/2104 to 21.29−0.02τ− 33.79e −τ/254.3 .This pattern is repeated at other sites, with the exponential timescales tending to be shorter when the simultaneous fit is done.For most deployments the timescale is less than 60 days.
This simultaneous fit to all the deployments gives an annual mass cycle m o with amplitude 0.86 mbar, phase 277 • , peaking at 8 October.

Noise on the simultaneous fit
We test the sensitivity to noise of the simultaneous fittings.We do this with a separate noise signal added to every deployment.The resulting distribution is shown in Fig. 8b and again as the ECCO model in Fig. 12.The annual cycle of m o has 95 % of results within amplitude range 0.68-1.05mbar, phase range 265-290 • (maximum at 26 September-21 October).
Figure 10 shows the comparison between the remainder p rec -p drift -p dyn -p t -p a ,where the ocean dynamics are from ECCO, and F s (m o ), where F s is taken from GLDAS-1 and m o is the peak of this distribution.We can see that at most sites most of the monthly or longer variability is accounted for.The exceptions are site 2 (Bay of Bengal), for which there is interannual variability with a larger signal in 2008, and sites 1 (Caribbean Sea) and 7 (Peru Basin), for which there is unaccounted variability of several mbar over periods of a few months.This variability is not captured by our noise model, effectively giving too high weighting to these sites.

Sensitivity to data selection
As a third method of estimating errors, we use bootstrappingwith-replacement.We select half of the deployments at random, and fit the annual to these simultaneously as described in Sect.5.1.This is repeated 100 times.The amplitude and phase of the annual is shown in Fig. 11.Note that bootstrapping should overestimate the errors as we are not using all the available deployments, and the sampling will sometimes select short deployments.
Also in Fig. 11 is shown the result for the simultaneous fit to only sites 11 and 17 (S and N from Hughes et al., 2012) and excluding those sites.The change between the latest calculation from S and N only (blue up-triangle) and the result quoted by Hughes et al. (2012) is largely due to the inclusion of the SAL corrections to the dynamical ocean pressure.The more detailed long-period tides also make a slight difference.All these results are enclosed by the 95 % noise contour for sites 11 and 17 only.The annual cycle of m o for sites 11 and 17 only has an amplitude 0.81 mbar, phase 272 • peaking 2 October.For all other sites: amplitude 0.89 mbar, phase 280 • peaking 10 October.Although the most probable value is not much changed from just using two sites, the error bounds are reduced by including more deployments.

Minimum deployment length
It would seem likely that it is easier to distinguish between drifts and annual cycles in long deployments.To test this we apply the simultaneous dedrifting procedure described in Sect.5.1 to bottom-pressure records, omitting deployments shorter than 6, 12, 18 or 24 months.We find that there is little difference in the amplitude or phase of the fitted annual ocean mass.This is perhaps because the few long deployments contain the majority of the data points, dominating the fit.The annual cycle m o found for deployments of minimum length [6,12,18,24] months has amplitude [0.86, 0.87, 0.88, 0.88] mbar and phase [279,279,282,277] degrees, peaks at [09,10,13,08] October.
We also tried applying the dedrifting procedure to individual sites omitting short deployments.Again, the long deployments seem to dominate, and the only deployments we have at > 2 years, at sites 2, 6, 10 and 12, all give similar results to the full records at those sites.
If only deployments shorter than 14 months are used, there is no convergence to an annual cycle, but the simultaneous fit to multiple sites will converge with only deployments shorter than 16 months (to amplitude 0.87 mbar, phase 270 • , peaking 30 September).

Comparison of hydrology models
If we use GLDAS-2.0 instead of GLDAS-1, and fit simultaneously to all sites (with ECCO), the ocean mass has an annual amplitude of 0.83 mbar, peaking at 279 • (9 October).This is a slightly smaller amplitude than with GLDAS-1, but the change is much less than when we used the technique of combining fits to individual sites.We think that the method of fitting to all sites simultaneously is less sensitive to the change in F s at site 2 between hydrology models.

Comparison of ocean models
Figure 12 shows the result of fitting to all sites simultaneously, using the various models available.The error bars are much larger for the OCCAM model, and NEMO4, which have much shorter overlap with the data (see Fig. 2).
Unlike the method of combining sites after fitting, with the simultaneous fitting the NEMO12 model increases the amplitude of the annual relative to the fit with the ECCO model.Again, the fit is less sensitive to site 2.
Summing the results of the ECCO and NEMO 1/12 models, which have the longest overlap with the data, gives a peak amplitude of 0.92 mbar and phase of 273 • (4 October), with 95 % of results within 0.71-1.14mbar or 262-288 • (22 September-19 October).

Selection of optimal sites
We also tried the calculation using only the "best" sites 6, 10-15 and 17.That is we excluded sites with less than 15 months records; sites 1 and 7, where there is unaccounted variability shown in Fig. 10; and site 2, where there are inconsistencies between GLDAS models and between ocean models.
For ECCO alone this moves the predicted annual to 0.82 mbar, peaking on 12 October.When NEMO12 is used there is very little change in the prediction, and if the models are combined then the annual has amplitude 0.85 mbar   and peaks at 10 October.The change is largely due to the omission of site 2.This estimate and probability distribution function is summarised and compared with results of other studies in Fig. 13.

Conclusions
We have shown that the annual cycle in ocean mass can be robustly determined from ocean bottom-pressure records of sufficient length.Records much shorter than a year contribute little to this determination, and for typical record lengths of 1-2 years, accurate calculation relies heavily on performing an iterative fit to the instrumental drift and annual cycle together.
The result from individual sites with records shorter than 16 months have unacceptably large error bounds, but by combining records from multiple locations in a simultaneous fit, we can provide a single estimate making using of all available data.The simultaneous fit is fairly robust to the selection of deployments, and even to omitting all deployments longer than 16 months.However, it does change slightly when only the "best" sites are used, omitting those with uncertainties between dynamical and GLDAS models, and we would recommend careful selection of sites.

Figure 2 .
Figure 2. Bottom-pressure data from the 17 test sites (mbar, offset for clarity), following subtraction by least squares fitting of tides with periods of fortnightly and shorter.Stars indicate the start of sensor deployments.There may be gaps in data within deployments.Black squares indicate convergence warnings in the drift fitting.These are not necessarily on the shortest deployments.The overlaid green line shows the sensor drift fitted using the iterative technique on all sites simultaneously.The end dates of the ocean models are also indicated.The y-axis offset between sites is 15 mbar.

Figure 3 .
Figure 3. Drift fitting to site 15.The yellow line is the fit to the raw bp record (blue), without taking account of the annual mass signal or other corrections.Green is the iterative fit to just that site, and magenta is the fit with the annual fitted to all sites simultaneously (discussed later).

Figure 4 .
Figure 4. Phasor diagrams for annual cycles at each site, for: (a) local bottom pressure p ANN m predicted using iterative fitting; (b) local bottom pressure p ANN m predicted without iterative fitting (sites 3 and 8 are outside the axes); (c) local pressure p h and global ocean mass m h predicted by two hydrology models: GLDAS-1 is shown with stars (p h ) and cross (m h ), GLDAS-2.0 with squares (p h ) and diamond (m h ); (d) global ocean mass m o , using a hydrology model GLDAS-1 to provide F s .Axes for (c) are indicated by red box on (a).All converged sites are included, some have large error bounds to be described below.

Figure 5 .
Figure 5. Scatter of p ANN m for each site with noise added to the bottom-pressure records.Subplot axes are identical to Fig. 4a.

Figure 6 .
Figure 6.Comparison of p ANN m derived from different models, for each site.Models are indicated as ECCO (blue), NEMO12 (green), NEMO4 (red), OCCAM12 (cyan), OCCAM4 (magenta).95 % of results lie within the contours.No contours are plotted if there is less than 13 months overlap between model and record data at a site.Subplot axes are identical to Fig. 4a.

Figure 7 .
Figure 7. Probability density functions (PDFs) for the amplitude and phase of m o calculated from each site with noise added to the bottompressure records, and mappings F s from hydrology and atmosphere model.Convergence was too poor at site 9 to calculate a PDF; 95 % of results fall inside the white contour.Subplot axes are identical to Fig. 4a.
results in the range of m o results shown in Fig. 4d.The most significant change is to site 2. With this conversion factor, we can convert the scattered estimates of p ANN m shown in Fig. 4a into corresponding estimates of m o in Fig. 4d.

Figure 8 .
Figure 8.(a) The probability of the annual mass having a given amplitude and phase for all 17 sites, combining individual site fits, renormalised to integrate to 1.(b) Scatter of annual amplitude and phase for all sites simultaneously, with noise added to the bottom-pressure records.Both plots are for the ECCO model, with GLDAS-1.95 % of results fall inside the white contour.

Figure 9 .
Figure 9. Amplitude of F s = p h /m h as derived from GLDAS-1 and GLDAS-2.0 for each site.

Figure 10 .
Figure10.Monthly mean bottom pressure at each site after removing ocean dynamics (from ECCO), tides, atmospheric pressure; vs. siteadjusted annual mass p ANN m calculated from simultaneous fit to bottom-pressure recorder data using ECCO with GLDAS-1 (peak of noise distribution).Bottom pressure before averaging is plotted in grey at each site.The y-axis offset between sites is 3 mbar.

Figure 11 .
Figure 11.Spread of simultaneous results under bootstrapping (dots).Also the result for only sites 11 and 17 (S and N from previous study) and excluding sites 11 and 17; the 95 % noise contour for only sites 11 and 17; the result from Hughes et al. (2012); and only sites 11 and 17 with no self-attraction and loading correction to the ocean dynamics.

Figure 12 .
Figure12.95 % noise contours for annual amplitude and phase fitted to each site simultaneously with noise, using ECCO (blue), NEMO12 (green), NEMO4 (red), OCCAM12 (cyan) or OCCAM4 (magenta).Note that the latter models have only a short overlap in time with the bottom-pressure data.

Figure 13 .
Figure 13.Probability distribution function of annual amplitude and phase for best sites (6, 10-15, 17) simultaneously with noise added to the bottom-pressure records, using both ECCO and NEMO12.95 % of results fall inside the white contour.The peak of the distribution is emphasised with a black cross, and grey dashed lines indicate one standard error and the 95 % bounds for amplitude and phase independently.With predictions for the ocean mass annual from previous studies: colours indicate publication, basis of method is indicated by circles (GRACE), squares (altimetry-steric), diamonds (ECCO model), stars (hydrology), up-triangles (bottom pressure), down-triangles (GRACE with complementary constraints including bottom-pressure data).Note that axes are zoomed from earlier figures, PDF colours are the same.

Table 1 .
Ocean circulation models used in this study.