Water transport among the world ocean basins within the water cycle

Global water cycle involves water-mass transport on land, atmosphere, ocean, and among them. Quantification of such transport, and especially its time evolution, is essential to identify footprints of the climate change and helps to constrain and improve climatic models. In the ocean, net water-mass transport among the ocean basins is a key, but poorly estimated parameter presently. We propose a new methodology that incorporates the time-variable gravity observations from 10 the GRACE satellite (2003-2016) to estimate the change of water content, and that overcomes some fundamental limitations of existing approaches. We show that the Pacific and Arctic Oceans receive an average of 1916 (95% confidence interval [1812, 2021]) Gt/month (~0.72 ± 0.02 Sv) of excess freshwater from the atmosphere and the continents that gets discharged into the Atlantic and Indian Oceans, where net evaporation minus precipitation returns the water to complete the cycle. This is in contrast to previous GRACE-based studies, where the notion of a seesaw mass exchange between the Pacific and 15 Atlantic/Indian Oceans has been reported. Seasonal climatology as well as the interannual variability of water-mass transport are also reported.

and Pacific Oceans (Warren, 1981). The Indian Ocean returns saltier water, but Pacific and Arctic Oceans return less-salty 30 waters, producing a salinity imbalance in the Atlantic. To restore the balance, freshwater must be transported outside the Atlantic at the rate of 0.13-0.32 Sv through the atmosphere (Zaucker et al., 1994). This WT produces an excess of freshwater in other ocean regions, as in the Pacific and Arctic Oceans, that must discharge out through the ocean.
Meanwhile, conventional observations on the lateral WT of world ocean climatology have been sparse. In fact, measuring such WT in an open ocean region proves difficult as it amounts only to a few tenths Sv, several orders of 35 magnitude smaller than the total ocean water inflow/outflows in such regions. For example, the Pacific is believed to receive regularly an inflow of 157 ± 10 Sv to south of Australia (Ganachaud & Wunsch, 2000), against three outflows: 0.7-1.1 Sv through the Bering Strait (Woodgate et al., 2012), 16 ± 5 Sv through the Indonesian Strait (Ganachaud & Wunsch, 2000), and 140-175 Sv through the Drake Passage (Ganachaud & Wunsch, 2000;Donohue et al., 2016;Colin de Verdière & Ollitrault, 2016;Vigo et al., 2018). 40 In this work we propose a new methodology devised to estimate the net WT through the boundaries of a given oceanic region. A defining feature of the proposed approach is the use of the time-variable gravity data from the GRACE (Gravity Recovery and Climate Experiment) satellite mission to estimate de change of water content. We apply the methodology, in conjunction with conventional meteorological data of general hydrologic budget schemes, to estimate the 45 time evolution over the period 2003-2016 of the net WT and exchanges among the four major ocean basins -namely Pacific, Atlantic, Indian, and Arctic. We analyse and report our results of the seasonal climatology as well as the interannual variability of WT. Such information, not available previously, would help elucidate the role of the oceans within the water cycle, and constrain ocean models (Warren, 1983;Rahmstorf, 1996;Emile-Geay et al., 2003;de Vries & Weber, 2005;Dijkstra, 2007). 50

Methodology and Data
The general hydrologic budget equation states that, at any continental location and any moment in time, the change of water content dW equals the precipitation P minus evapotranspiration E (as vertical transport) minus the net runoff R (as horizontal transport): Under the conservation of water mass, the global net P−E over ocean is negative [e.g., Hartmann, 1994]. That amount of water gets transported to land through atmosphere and returns to the ocean as R completing the water cycle. The general R for a river basin connected to the ocean consists of river runoff, land ice melting, and submarine groundwater discharge to 60 ocean.
For an ocean region, R represents the inflow from adjacent land regions plus an extra additive term, call it N, accounting for water exchange between neighbouring ocean regions through boundaries, as (positive) inflow or (negative) outflow: 65 The ocean water flux N is the target quantity that we shall solve for as a residual in Equation 2, which up till now has been infeasible to estimate directly [Rodell et al., 2015]. Note that N represents the integrated WT over the total-column depth of 70 ocean, including deep-water flows. This is a strength of the GRACE observation for the oceans, compared to in-situ or other remote-sensing measurements typically targeting only the surface layer.
Our targeted four ocean basins are largely separated geographically with designated continental boundaries and restricted water throughways. The land is divided into their associated drainages according to the global continental runoff 75 pathways scheme of Oki & Sud (1998). There are no direct water exchanges in the form of R among land drainages (see Figure 1). The WT component R is estimated through Equation (1) over each continental region, then input to Equation (2) to estimate N in the associated ocean basins.
The P-E data we use are adopted from the ERA5 reanalysis [Hersbach, 2018], which assimilates observations into 80 general-circulation modelling provided by the European Centre for Medium-Range Weather Forecasts (ECMWF). They are given at 0.25º latitude/longitude regular grids and monthly (and hourly) intervals for global coverage of both continents and oceans. In order to match the spatial resolution of the above-mentioned continental runoff pathways data, we homogenise the grid to 1ºx1º by averaging the corresponding 0.25º grid points.

85
The critical knowledge needed in Equations (1) and (2), now obtainable from GRACE monthly data, is dW (Tapley et al., 2004(Tapley et al., , 2019, the month-to-month difference of the stored water. Note that the GRACE mass variability pertains to WT directly, as opposed to, for example, altimetric sea level measurements that also contain non-WT, steric effects. We use the GRACE "mascon" (mass concentration) solutions that have already been converted into surficial mass from the original time-variable gravity observations (in our case the GRACE RL06 mascon dataset provided by the Center of Space Research 90 of University of Texas; see Save et al., 2016, Save, 2019. The non-surficial gravity change due to the glacial isostatic adjustment (GIA) has been removed to the extent of the ICE6G-D model (Peltier et al., 2018). Any other non-surficial effect such as long-term tectonics do masquerade as mascons (Chao, 2016) but are here ignored; so are the non-climatic sources such as the rare, local earthquake events. As is a common practice, the C20 Stokes coefficient is replaced with the solution from Satellite Laser Ranging (Cheng & Ries, 2017). GRACE is not sensitive to the geocenter variations, and its degree-1 Stokes coefficients are set to zero. We had tried adding to GRACE data an estimate of geocenter variations due to modelled water-mass variability (Swenson et al., 2008), and our reported results would change less than 1%. On the other hand, we do add back the GAD product, the oceanic effects on gravity change according to the operational numerical weather prediction (NWP) model from ECMWF and to an unconstrained simulation with the global ocean general circulation model MPIOM, respectively (Dobslaw et al., 2017). The GAD had beforehand been removed from the processing of the GRACE data for de-100 aliasing purposes, now added back so that we recover the "true" ocean mass variability. Data are provided on a 0.25º regular grid; we reduce it to 1º regular grids, still finer than the spatial resolution of GRACE (~300 km), to match the spatial resolution of the continental drainage basin data as above.
GRACE's degree-0 Stokes coefficients DC00 is set identically to zero on the recognition that the total mass is 105 constant. The GAD product adjustment is performed by the GRACE project in such a way that the atmospheric and dynamic oceanic mass changes are removed beforehand. After that, the GRACE DC00 are still set to zero even though they should match the opposite of the removed signals by GAD. To restore the lost degree-0 signal, the GAD product must be added back to GRACE with averaged ocean signal set to zero, and then, the DC00 from an atmospheric model must be subtracted from GRACE data. Doing so, the GRACE data will account for the global exchange of water-mass between the Earth 110 surface and atmosphere. Such correction has recently proved to improve the agreement between the GRACE global ocean mass change and non-steric sea level variation estimates from altimetry and ARGO data (Chen et al., 2019). Looking for consistency between the GRACE and ERA5 datasets, we use DC00 from P -E to restore degree-0 signal in dW. This DC00 accounts for uniform mass variations in the global surface equivalent to a global averaged signal for P−E, at 188 Gt/month (95% confidence interval CI95=[136, 243], see below). As global -(P−E) represents the variability of global total-column 115 water (TCW), it should match the time derivative of the global TCW. However, the average rate of change of the global TWC in ERA5 is 1.5 Gt/month (CI95=[−9.2, 12.7]), although in the range of previously reported values of [−0.9, 4.3] Gt/month [Nilsson & Elgered, 2008] departs far from the global -(P−E) value. This reveals some internal inconsistency within the ERA5 dataset. However, while artificially increasing the dW estimate, the excessive value of P−E does not affect the WT components R and N estimated from Equations (1) and (2), since the degree-0 signal vanishes due to the residual 120 estimate between dW and P−E. In fact, adding DC00 from P−E to GRACE is numerically equivalent to setting P−E DC00 to zero as far as Equations (1) and (2) are concerned.
The reported 95% confidence intervals and the correlation coefficients are evaluated using the stationary bootstrap scheme of Politis and Romano (1994) (with optimal block length selected according to Patton et al., 2009), and the quantile 125 method. The number of bootstrap replications is set to 2000. In general, half length of the confidence interval can be very well approximated by twice the standard deviation of the sample mean estimated from the bootstrap replications. Prior to applying the bootstrap, least-squares estimated linear/quadratic trend and sinusoid with the most relevant frequencies are https://doi.org/10.5194/esd-2020-54 Preprint. Discussion started: 20 July 2020 c Author(s) 2020. CC BY 4.0 License. removed to meet the stationarity conditions of the method and to avoid spurious correlation. In particular, for the WT N component we proceed as follows: (i) a model with linear, annual, and semiannual signals is fitted to the data and subtracted; 130 (ii) 2000 bootstrap samples of the residuals are generated; (iii) these are added back to the fitted model to obtain 2000 bootstrap samples of the original time series. These samples are then used to obtain point estimates of the mean fluxes, the annual amplitudes and phases, and the corresponding quantile-based 95% confidence intervals. In order to study the robustness of the results with respect to the model choice, the analysis is rerun using 11 alternative models obtained considering different forms for the trend component (quadratic or constant) and including higher frequencies in the harmonic 135 regression (up to 5). The results are robust. The relative difference with respect to the reported values is smaller than 1.2% for point estimates and smaller than 3.3% for the extremes of the 95% confidence intervals.

Results
The various WT components of the Pacific and its associated land drainage regions are shown in Figure 2 in units of 140 Gt/month (1 Sv ≈ 2600 Gt/month, 1 Gt = 10 12 kg). The same analysis is applied to the rest of the ocean basins, i.e. the AIA oceans individually and collectively, with its associated land drainages (see Figure 1).

Mean fluxes
Averaged over the studied 14 years, the Pacific loses water through the atmospheric P−E at the average rate of 142 145 Gt/month (CI95=[48,243]), which is greatly over-compensated by inflow R from land of 1403 Gt/month (CI95= [1370,1436]).

150
In the AIA Oceans, the situation is found to be markedly distinct, given the fact that the AIA oceans together have surface area comparable to the Pacific (177x10 6 m 2 which will be called the "AIA inflow". As expected from the overall conservation of water mass inherent in our methodology, the estimated Pacific outflow and AIA inflow match (Figure 3). It is worth mentioning that a difference of 188 Gt/month would exist between the two mean flux values if the degree-0 correction were not applied.

Annual climatology 175
The WT climatology of the N component is estimated in two ways: (1) averaging the 14 N values for each months of the year (Figure 5a); and (2) fitting a linear trend plus annual and semiannual components model as described in Section 2. Annual amplitudes and phases are plotted in Figure 5b and reported, with corresponding 95% quantile-based confidence intervals, in Table 1.

180
The Pacific and Arctic Oceans show an overall outflow throughout the year, unlike the Atlantic and Indian Oceans, which show an inflow for every month. The Pacific outflow shows a prominent seasonal undulation peaked around August 3 and a peak-to-peak WT variation of ~2000 Gt/month from boreal summer to November, when a near-zero minimum occurs. The Arctic Ocean show half of the Pacific variability and a less pronounced seasonal undulation. A minimum outflow of ~320 Gt/month is reached in March and April, and a maximum ~1300 Gt/month in July. Together, the Pacific and Arctic Oceans 185 send ~3000 Gt/month of seawater to the Atlantic and Indian Oceans during boreal summer, and a minimum amount five times lower, around 600 Gt/month, in November. The annual maximum is reached on August 8th. The Atlantic/Arctic inflow shows a specular behaviour. Separately, the Atlantic and Indian inflows show a similar peak-to-peak variation of ~2000 Gt/month, reaching the maxima in August and May, respectively. The Indian maximum seems to be related to a local maximum of the Pacific outflow. The annual maxima of net WT of the four basins are reached between August 3 rd and 190 September 9 th , although the annual signals of the Pacific and Indian Oceans almost triple those from Arctic and Atlantic Oceans (Table 1 and Figure 5b). https://doi.org/10.5194/esd-2020-54 Preprint. Discussion started: 20 July 2020 c Author(s) 2020. CC BY 4.0 License.

Interannual variability
Interannually, the Pacific outflow shows remarkable variability, mainly produced by P on the continents, which is 195 inherited by R, and P−E in the oceans (Figure 2 Figure 6 show all indices with amplitudes normalized to one standard deviation, as well as the de-trend, de-season, standard deviation normalized Pacific outflow. The correlation analysis between the Pacific 215 outflow and the SOI shows no overall correlation during the whole period, meaning that the influence of ENSO on the Pacific outflow may be restricted to the strong phases of ENSO as in 2009 and 2010. A similar lack of correlation is observed with respect to the AMO, AAO, and AO.

Discussion and Conclusions. 220
In this work we present a new methodology that combines GRACE data with the general hydrologic budget equation to estimate the horizontal water-mass convergence/divergence for any oceanic region. We have assumed that the gravity changes are produced by mass changes on the Earth surface, such as in the oceans, so that the mascon solution is physically meaningful (Chao, 2016). Any mis-modelling of the ocean basin "container" volume change due to GIA and other non-https://doi.org/10.5194/esd-2020-54 Preprint. Discussion started: 20 July 2020 c Author(s) 2020. CC BY 4.0 License.
surficial changes would masquerade as WT variations. However, they are not critical as far as our non-secular analysis is 225 concerned.
We use the proposed methodology to estimate the net WT and exchanges among the Pacific, Atlantic, Indian, and Arctic Oceans, for the period of 2003 -2016. Our main finding is that the Pacific and Arctic Oceans, while replenished with precipitation and land runoff, are nearly continuously losing water to the Atlantic and Indian Oceans. In particular, the WT 230 climatology is such that the Pacific Ocean loses water at a rate from near zero to up to the peak of 2000 Gt/month during the boreal summer, which coincides with the maximum of the global atmosphere water content. On top of the climatology, the interannual Pacific water loss varies significantly between ~950 to ~1450 Gt/month annual means during the studied period, but seemingly uncorrelated with ENSO.

235
The results presented here are consistent with the well-known salinity asymmetry between the Pacific and Atlantic Oceans (Reid, 1953;Warren, 1983;Broecker et al., 1985;Zaucker et al., 1994;Rahmstorf, 1996;Emile-Geay et al., 2003;Lagerloef et al., 2008;Czaja, 2009;Reul, 2014). However, they are in contrast to previous GRACE-based studies where a simple seesaw WT between the Pacific and the Atlantic/Indian oceans was reported [Chambers & Willis, 2009;Wouters et al., 2014]. In those studies, the P−E+R term in Equation 2 in both Pacific and Atlantic/Indian Oceans was approximated by that 240 from the global ocean mean. However, the mean freshwater flux in the Pacific (1261 Gt/month) quite mis-matches that in the Atlantic/Indian Oceans (−1837 Gt/month), meaning that the approximation was quite poor and hence the N term was not properly estimated in these studies (see Supplementary Materials for further discussion).
Differences in freshwater fluxes between the Pacific and Atlantic Oceans produce salinity contrasts, and in turn contrasts on 245 deep water formation. Nevertheless, there are other factors influencing these contrasts such as the narrower extent of the Atlantic (de Boer et al., 2008;Jones & Cessi, 2017), the meridional span of the African and American continents (Nilsson et al., 2013;Cessi & Jones, 2017), and the salty WT from the Indian Ocean to the Atlantic (Gordon, 1986;Marsh et al., 2007). AMOC is also influenced by WT through Bering Strait (Reason & Power, 1994;Goosse et al., 1997;Wadley & Bigg, 2002), and by surface processes of temperature, precipitation and evaporation at low-latitudes of Pacific and Indian Oceans 250 (Newsom & Thompson, 2018). The relative importance among the multiple drivers influencing the AMOC is an open problem (Ferreira, 2018). The net WT estimated here provides information for differences between oceanic inflows and outflows, which can be useful to elucidate on this problem.
Net WT in the open oceans can alternatively be estimated using global ocean models, which simulates observational data 255 based on physical principles. However, these models are not necessarily sensitive to the WT specifically given the data types, and the geography and topography resolutions involved in the models. Knowledge about three-dimensional global ocean circulation could also elucidate on the net WT. However, the small ratio between the net and the total WT hinders the estimation of the former from the latter.

260
We have applied our WT estimation scheme to the four major ocean basins. The methodology can of course be applied to any extensive ocean region of interest as long as it is much larger than the GRACE resolution. The findings reported here will be useful for a better understanding of the global climate system in terms of its climatology and spatio-temporal variations.

Appendix: Apparent net mass exchange between Pacific and Atlantic/Indian oceans
We shall here that the net water mass exchange between the Pacific and Atlantic/Indian Oceans reported by Chambers and Willis (2009) was a mathematical artefact. Their Equation (2)  hence their resultant estimates of WT are rather poor. In addition, under their approximation an apparent net mass exchange will always arise, since Thus, wherever the signal is in the Pacific and Atlantic/Indian oceans, the anomalies with respect to the global ocean mean will always mirror each other, showing an apparent net mass exchange between them, even if such exchange does not exist.   https://doi.org/10.5194/esd-2020-54 Preprint. Discussion started: 20 July 2020 c Author(s) 2020. CC BY 4.0 License.