Water transport among the world ocean basins within the water cycle

The global water cycle involves water-mass transport on land, in the atmosphere, in the ocean, and among them. Quantification of such transport, especially its time evolution, is essential to identify the footprints of climate change, and it also helps to constrain and improve climatic models. In the ocean, net water-mass transport among the ocean basins is a key process, but it is currently a poorly estimated parameter. We propose a new methodology that incorporates the time-variable gravity observations from the Gravity Recovery and Climate Experiment (GRACE) satellite (2003–2016) to estimate the change in water content; this new approach also overcomes some fundamental limitations of existing methods. We show that the Pacific and Arctic oceans receive an average of 1916 (95 % confidence interval of [1812, 2021]) Gt per month ($\sim \mathrm{0.72}\pm \mathrm{0.02}$ Sv) of excess freshwater from the atmosphere and the continents that is 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 see-saw mass exchange between the Pacific and the Atlantic and Indian oceans has been reported. Seasonal climatology as well as the interannual variability of water-mass transport are also reported.

1 Introduction

The water-mass transport (hereafter WT for brevity) in the oceans is a deciding factor in the world climate system. Quantification of such transport, especially its time evolution, is essential to better understand climate change. The Atlantic Ocean presents a notable freshwater flux deficit, in contrast to the Pacific Ocean. This produces a salinity asymmetry that explains why deep waters are formed in the North Atlantic and not in the North Pacific (Warren, 1983; Broecker et al., 1985; Rahmstorf, 1996; Emile-Geay et al., 2003; Czaja, 2009). The upper layers of the North Atlantic flow northward, while deep waters flow southward, forming the Atlantic Meridional Overturning Circulation (AMOC), which distributes heat within the Earth system and influences temperature and precipitation patterns worldwide (Vellinga and Wood, 2002). While small changes in the hydrological cycle may have caused changes in AMOC during the last glaciation that led to abrupt climate changes (Clark et al., 2002), different models project a weakening of the AMOC in the 21st century that would lead to profound climatic and ecological changes over vast regions (Collins et al., 2019). The Antarctic Circumpolar Current (ACC) receives deep water injected by the AMOC with excess salinity, which in turn is transported into the Indian and Pacific oceans (Warren, 1981). The Indian Ocean returns saltier water, but the Pacific and Arctic oceans return less-salty waters, producing a salinity imbalance in the Atlantic. To restore the balance, freshwater must be transported outside the Atlantic at a 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 of a sverdrup (Sv), which is several orders of magnitude smaller than the total ocean water inflows and outflows in such regions. For example, the Pacific is believed to regularly receive an inflow of 157±10 Sv to the south of Australia (Ganachaud and Wunsch, 2000), against three outflows: 0.7–1.1 Sv through the Bering Strait (Woodgate et al., 2012), 16±5 Sv through the Indonesian straits (Ganachaud and Wunsch, 2000), and 140–175 Sv through the Drake Passage (Ganachaud and Wunsch, 2000; Donohue et al., 2016; Colin de Verdière and Ollitrault, 2016; Vigo et al., 2018).

In this work, we propose a new methodology that has been 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 the change in water content. We apply the methodology, in conjunction with conventional meteorological data of general hydrologic budget schemes, to estimate the time evolution of the net WT and exchanges among the four major ocean basins – namely the Pacific, Atlantic, Indian, and Arctic – over the period from 2003 to 2016. We analyse and report our results of the seasonal climatology as well as the interannual variability of WT. Such information, which has not been previously available, is of valuable importance. For example, in closed regions, net WT through the boundaries on the surface must be counteracted by moisture fluxes through the same boundaries in the atmosphere in order to match GRACE measurements. Such an approach has been successfully applied to study the hydrological cycle of South America (Liu et al., 2006). At the ocean basin scale, knowledge about net WT would not only help elucidate the role of the oceans within the water cycle, but it would also impose restrictions on moisture advection in the atmosphere that would help to improve atmospheric models. Moreover, ocean models usually deal with inflows and outflows of a given ocean region (Warren, 1983; Rahmstorf, 1996; Emile-Geay et al., 2003; de Vries and Weber, 2005; Dijkstra, 2007), and net WT estimates for such ocean regions would be useful to impose constraints on the relationship between their inflows and outflows, which would improve the reliability of the models. Better models will improve our knowledge of the Earth's WT dynamics and its evolution in the future, which is critical in the present scenario of climate change.

2 Methodology and data

2.1 Methodology

The general hydrologic budget equation states that, at any continental location and any moment in time, the change in water content (dW) equals the precipitation (P) minus evapotranspiration (E; as vertical transport) minus the net runoff (R; as horizontal transport):

$$\begin{array}{}\text{(1)}& dW=P-E-R,\mathrm{for}\mathrm{land}.\end{array}$$

Under the conservation of water mass, the global net P−E over ocean is negative (e.g. Hartmann, 1994). That amount of water is transported to land via the atmosphere and returns to the ocean as R, thereby completing the water cycle. The general R for a river basin connected to the ocean consists of river runoff, land ice melt, and submarine groundwater discharge to ocean. The R component will be estimated as a residual in Eq. (1).

For an ocean region, R represents the inflow from adjacent land regions as well as an extra additive term, N, which accounts for water exchange between neighbouring ocean regions through boundaries as (positive) inflow or (negative) outflow:

$$\begin{array}{}\text{(2)}& dW=P-E+R+N,\mathrm{for}\mathrm{ocean}.\end{array}$$

The ocean water flux N is the target quantity that we shall solve for as a residual in Eq. (2), which until now has been infeasible to estimate directly (Rodell et al., 2015). Note that N represents the integrated WT over the total-column depth of ocean, including deep-water flows. This is a strength of the GRACE observation for the oceans, compared with in situ or other remote-sensing measurements that typically only target 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 pathways scheme of Oki and Sud (1998). There are no direct water exchanges in the form of R among land drainages (see Fig. 1). The WT component R is estimated via Eq. (1) over each continental region and then input to Eq. (2) to estimate N in the associated ocean basins.

Figure 1 The Pacific, Atlantic, Indian, and Arctic ocean basins and their associated continental drainage basins according to the global continental runoff pathways scheme of Oki and Sud (1998). Within each basin, the darker colour represents the continental basin, and the lighter colour represents the ocean basin. White regions represent endorheic basins.

2.2 Precipitation and evaporation data

The P and E data that we use are adopted from the ERA5 reanalysis (Hersbach et al., 2018), which assimilates observations into 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 homogenize the grid to $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ by averaging the corresponding 0.25^∘ grid points.

2.3 Time-variable GRACE data

The critical knowledge needed in Eqs. (1) and (2), which is now obtainable from GRACE monthly data, is dW (Tapley et al., 2004, 2019) – the month-to-month difference in the stored water. Note that the GRACE mass variability pertains directly to WT, 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, CSR, 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 would be incorrectly interpreted as water-mass fluxes (Chao, 2016), but they may only be important in the determination of secular trends. Additional non-climatic sources, such as the rare local earthquake events, may be similarly misinterpreted. As the C₂₀ Stokes coefficient is not well determined from the GRACE mission, it is replaced by a more accurate solution from satellite laser ranging (SLR; Cheng and Ries, 2017). GRACE is not sensitive to geocenter variations, and its degree-1 Stokes coefficients are set to zero. We tried adding an estimate of geocenter variations due to modelled water-mass variability to GRACE data (Swenson et al., 2008) and found that our reported results only changed by less than 1 %. Furthermore, the atmospheric effects as well as some oceanic effects on gravity change have already been removed from the processing of the GRACE data by CSR, for de-aliasing purposes, according to the operational numerical weather prediction (NWP) model from ECMWF and to an unconstrained simulation with the global ocean general circulation model MPI-OM (Max-Planck-Institute Global Ocean/Sea-Ice Model; Dobslaw et al., 2017). To recover the “true” ocean mass variability, we restore the removed signal on the oceans by adding back the GAD product, which is set to zero on the continents. Data are provided on a 0.25^∘ regular grid; we reduce the resolution to 1^∘ regular grids, which is still finer than the spatial resolution of GRACE (∼300 km), in order to match the spatial resolution of the continental drainage basin data as above.

GRACE's degree-0 Stokes coefficients ΔC₀₀ are set to be identical to zero on the basis that Earth's total mass (including the atmosphere) is constant. Any increase (decrease) in the water-mass budget of the atmosphere will then be counteracted by a decrease (increase) of the same amount of water mass at the surface. However, after the atmospheric and dynamic oceanic mass changes are corrected for in GRACE data, the GRACE ΔC₀₀ values are still set to zero, even though they should match the opposite of the removed signals. To restore the lost degree-0 signal, the GAD product (which is set to zero on the continents) must be added back to GRACE with the averaged ocean signal set to zero; the ΔC₀₀ from an atmospheric model must then be subtracted from GRACE data in order to force the Earth's total mass to be constant. After these steps have been undertaken, the GRACE data will account for the global exchange of water mass between the Earth surface and the atmosphere. Such a correction has recently proven 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 ΔC₀₀ from P−E to restore the degree-0 signal in dW. This ΔC₀₀ accounts for uniform mass variations in the global surface that are equivalent to a global averaged signal for P−E, at 188 Gt per month (95 % confidence interval CI₉₅=[136, 243], see below). As global $-(P-E)$ represents the variability of global total-column water (TCW), it should match the time derivative of the global TCW. However, the average rate of change in the global TCW in ERA5 is 1.5 Gt per month (CI${}_{\mathrm{95}}=[-\mathrm{9.2}$, 12.7]); although this is in the range of previously reported values of [−0.9, 4.3] Gt per month (Nilsson and Elgered, 2008), it is far removed from the global $-(P-E)$ value. This reveals some internal inconsistency within the ERA5 dataset. However, while this artificially increases the dW estimate, the excessive value of P−E does not affect the WT components R and N estimated from Eqs. (1) and (2), as the degree-0 signal vanishes due to the residual estimate between dW and P−E. In fact, adding ΔC₀₀ from P−E to GRACE is numerically equivalent to setting P−EΔC₀₀ to zero as far as Eqs. (1) and (2) are concerned.

2.4 Confidence intervals

The reported 95 % confidence intervals and the correlation coefficients are evaluated using the stationary bootstrap scheme of Politis and Romano (1994) (with the optimal block length selected according to Patton et al., 2009) as well as the percentile method. The intuition underlying the bootstrap scheme is simple. Suppose that the observed time series x₁, …, x_n is a realization of the random vector (X₁, …, X_n) with the joint distribution P_n and which is assumed to be part of a stationary stochastic process. Given X_n=(X₁, …, X_n), we first build and estimate ${\widehat{P}}_n$ of P_n. B random vectors ($X_{\mathrm{1}}^{*}$, …, $X_n^{*}$) are then generated from ${\widehat{P}}_n$. If ${\widehat{P}}_n$ is a good approximation of P_n, the relation between ($X_{\mathrm{1}}^{*}$, …, $X_n^{*}$) and ${\widehat{P}}_n$ should reproduce the relation between (X₁, …, X_n) and P_n well (for an introduction to bootstrap methods for time series, see Kreiss and Lahiri, 2012, and the references therein). Here, the number of bootstrap replications was set to B=2000. In general, the half-length of the confidence interval can be approximated very well by twice the standard deviation of the sample mean estimated from the bootstrap replications. Prior to applying the bootstrap to a time series, the least-squares estimated linear and quadratic trend as well as the sinusoid with the most relevant frequencies are removed from the series to meet the stationarity conditions of the method. In particular, each series has been decomposed into trend, seasonal, and residual components. The bootstrap is applied to the residual component and produces bootstrap samples of the residuals. For the evaluation of confidence intervals for the different components of WT, the trend and seasonal terms are added back (to the bootstrap sample of the residuals), producing bootstrapped time series of the component of interest. These samples are then used for further analysis. As an illustration, for the WT N component, we proceed as follows: (i) a model with linear, annual, and semiannual signals is fitted to the data, and the fitted linear trend and the annual and semiannual signals are subtracted from the original time series; (ii) the stationary bootstrap is then applied to the residuals, producing 2000 bootstrap samples of the residuals; (iii) the estimated trend and seasonal components are added back to each bootstrap sample of the residuals which results in an ensemble of 2000 bootstrapped time series for the N component; (iv) these 2000 bootstrapped time series are used to obtain the 95 % confidence intervals for the mean fluxes (average of N over the 14-year period of study) and for the amplitude and phase of the annual component using the percentile method. For the mean fluxes, the average of N for each of the 2000 bootstrapped time series was first evaluated, and the 0.025 and 0.975 percentiles of these 2000 averages were then reported as the 95 % confidence interval. For the study of the climatology, a linear trend model with annual and semiannual components was fitted to the 2000 bootstrapped time series, producing corresponding estimates of the annual amplitude and phase. The 0.025 and 0.975 percentiles of these estimates were reported as the 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 that are obtained by considering different forms of the trend component (quadratic or constant) and including higher frequencies in the harmonic regression (up to 5). The results are robust. The relative difference with respect to the reported values is less than 1.2 % for point estimates and less than 3.3 % for the extremes of the 95 % confidence intervals.

As an independent check of the bootstrap, confidence intervals for the mean value of N have been also evaluated by propagating the error estimate in GRACE data (using the JPL GRACE mascon solution for which error estimates are available). The resulting intervals were consistent with those of the bootstrap method. In particular, we show that, in all cases, the bootstrap intervals contain the intervals obtained from error propagation (see Sect. 4 for details). In this respect, the CI₉₅ from bootstrap analysis can be considered a conservative estimate. This should be expected, as the residual component underlying the bootstrap approach includes measurement errors and other type of errors (related, for example, to the estimate of the trend and seasonal terms). As a result, the uncertainties in the transports estimated by the bootstrap should be larger than the corresponding uncertainties estimated by error propagation.

Note that the bootstrap was applied to the bivariate time series of the residuals of the two variables of interest for the study of correlation, producing an ensemble of 2000 bivariate time series of residuals. For each bivariate time series of residuals, the correlation between the two components of the series was first evaluated. The average and the 0.025 and 0.975 percentiles of these 2000 estimates were reported as a point estimate and confidence limits respectively for the correlation between the two variables of interest (the correlation between the residual components is used to avoid spurious correlation).

3 Results

The various WT components of the Pacific and its associated land drainage regions are shown in Fig. 2 in units of gigatonnes per month (1 Sv ≈ 2600 Gt per month; 1 Gt = 10¹² kg, the weight of 1 km³ of freshwater). The same analysis is applied to the rest of the ocean basins, i.e. to the “non-Pacific” AIA oceans, individually and collectively, with their associated land drainages (see Fig. 1).

Figure 2 WT from Eqs. (1) and (2) in the Pacific (a, b) and in the Atlantic, Indian, and Arctic (AIA) oceans collectively (c, d), as well as their drainage basins. (a, c) The associated land drainage basins and (b, d) ocean basins. Labels on the y axis correspond to the mean ± standard deviation of the associated curve. Thick lines are the low-pass-filtered signal by a Hann function of 24 months. All curves in the same panel are plotted on the same scale. P, E, and P−E are from the ERA5 dataset; dW is estimated from GRACE; R and N are estimated as a residual in Eqs. (1) and (2) respectively.

3.1 Mean fluxes

Averaged over the 14 years studied, the Pacific Ocean loses water through the atmospheric P−E at an average rate of 142 Gt per month (CI₉₅=[48, 243]), which is greatly overcompensated for by an inflow R from land of 1403 Gt per month (CI₉₅=[1370, 1436]). From this surplus, a minor (if any) amount of 67 Gt per month (CI₉₅=[25, 108]) remains (and accumulates) in the Pacific, while 1194 Gt per month (CI₉₅=[1096, 1291]) is transported horizontally to the “non-Pacific” Atlantic, Indian, and Arctic (AIA) oceans, which will be referred to as the “Pacific outflow” hereafter.

In the AIA oceans, the situation is found to be markedly distinct, given the fact that the AIA oceans have a combined surface area comparable to the Pacific (177×10⁶ km²). The AIA oceans collectively lose 3484 Gt per month (CI₉₅=[3406, 3560]) through the atmospheric P−E, which is ∼25 times the amount lost by the Pacific. Only ∼68 % of this water deficit is compensated for by a land R inflow of 2378 Gt per month (CI₉₅=[2337, 2419]). With the nominal minor amount of water accumulation of 87 Gt per month (CI₉₅=[44, 130]), the AIA oceans consequently present an average inflow of 1194 Gt per month (CI₉₅=[1102, 1284]) from the Pacific, which will be referred to as the “AIA inflow” hereafter.

As expected from the overall conservation of water mass inherent in our methodology, the estimated Pacific outflow and the AIA inflow match (Fig. 3). It is worth mentioning that a difference of 188 Gt per month would exist between the two mean flux values if the degree-0 correction was not applied.

Figure 3 Monthly time series of WT flux from the Pacific to the AIA oceans. The red curve is (the opposite of) the Pacific outflow, and the black curve is the AIA inflow. Thick lines are the low-pass-filtered signal by a Hann function of 24 months.

Figure 4 Same as Fig. 2 but for the Atlantic, Indian, and Arctic oceans.

Figure 5 Diagram of the mean values of the WT of the study regions. Units are in gigatonnes per month.

Corresponding analyses have been performed for the Atlantic, Indian, and Arctic oceans separately. The time evolution of the WT components in Eqs. (1) and (2) are shown in Fig. 4, and a diagram of the water-mass fluxes is shown in Fig. 5. On average, the Atlantic Ocean receives 926 Gt per month (CI₉₅=[876, 980]; or 0.36 Sv) of salty water, and it loses 879 Gt per month (CI₉₅=[828, 930]) to the atmosphere via $P-E+R$. The latter is equivalent to a freshwater deficit of 0.34 Sv, which increases the near-surface salt concentration and enables water to sink in the North Atlantic, producing deep water. These values are close to the 0.13–0.32 Sv estimated from ocean models, as needed to maintain the salinity balance in the Atlantic Ocean (Zaucker et al., 1994). Similarly, the Indian Ocean loses 957 Gt per month (CI₉₅=[894, 1022]) of freshwater that is restored by 991 Gt per month (CI₉₅=[907, 1073]) of salty water. The freshwater lost by the Atlantic and Indian oceans via $P-E+R$ goes to the Pacific (1261 Gt per month, CI₉₅=[1171, 1347]) and Arctic (730 Gt per month, CI₉₅=[712, 747]) oceans, which return 1194 (CI₉₅=[1096, 1291]) and 723 (CI₉₅=[708, 739]) Gt per month of salty water through the ocean respectively. Thus, the Pacific presents a surplus of freshwater that reduces the near-surface salt concentration, preventing the formation of deep water. Together, the Pacific and Arctic oceans supply 1917 Gt per month (CI₉₅=[1812, 2021]) of water to the Atlantic and Indian oceans, where it is reincorporated into the water cycle via net E−P. As in previous studies (see Craig et al., 2017, for a synthesis), the freshwater lost in the Indian Ocean is similar to that in the Atlantic Ocean. In these studies, $P-E+R$ is close to zero in the Pacific Ocean, producing a difference of 0.4 Sv between the Atlantic and Pacific oceans. In our study, $P-E+R$ is 1261 Gt per month in the Pacific Ocean, and the difference with the Atlantic increases to ∼0.8 Sv. Some of these differences would be expected, as the ocean basins are not defined in exactly the same way. On the other hand, the global R is 3781 Gt per month (or $\mathrm{3781}\times \mathrm{12}=\mathrm{45}\mathrm{368}$ km³ yr⁻¹), which is close to the 41 867 km³ yr⁻¹ reported by the Global Runoff Data Centre (GRDC, 2014). At the basin scale, R is 16 834 km³ yr⁻¹ in the Pacific, which is greater than the 11 826 km³ yr⁻¹ reported by GRDC. In the Atlantic, Indian, and Arctic, R is 18 228, 4479, and 5827 km³ yr⁻¹ respectively, which is closer to the respective GRDC values of 20 772, 5238, and 4080 km³ yr⁻¹. Finally, according to the diagram in Fig. 5, the water content in the atmosphere decreases by 178 Gt per month (and is gained by Earth's surface), but this amount is not realistic (as discussed in Sect. 2), as it should increase by a few gigatonnes per month (Nilsson and Elgered, 2008). This value differs from the 188 Gt per month mentioned in Sect. 2 because the endorheic regions are not considered here.

3.2 Annual climatology

The WT climatology of the N component is estimated in two ways: (1) averaging the 14 N values for each month of the year (Fig. 6a) and (2) fitting a linear trend as well as annual and semiannual components' model as described in Sect. 2. Annual amplitudes and phases are plotted in Fig. 6b and reported, with their corresponding 95 % quantile-based confidence intervals, in Table 1.

Figure 6 (a) Annual climatology time series (the error bar is 1 standard deviation) and (b) phasor diagram (amplitude in units of Gt per month, phase angle according to Eq. 3) of the inflow and outflow WT of the ocean basins.

Table 1 Mean and annual signals of the N component as estimated from the CSR mascon solution for different ocean basins according to Eq. (2).

The Pacific and Arctic oceans show an overall outflow throughout the year; this is in contrast to the Atlantic and Indian oceans, which show an inflow for every month. The Pacific outflow shows a prominent seasonal undulation that peaks around 3 August and a peak-to-peak WT variation of ∼2000 Gt per month from boreal summer to November, when a near-zero minimum occurs. The Arctic Ocean shows half of the Pacific variability and a less pronounced seasonal undulation. A minimum outflow of ∼320 Gt per month is reached in March and April, and a maximum outflow of ∼1300 Gt per month is seen in July. Together, the Pacific and Arctic oceans send ∼3000 Gt per month of seawater to the Atlantic and Indian oceans during boreal summer and a minimum amount 5 times lower (around 600 Gt per month) in November. The annual maximum is reached on 8 August. The combination of the Atlantic and Indian inflow mirrors this behaviour. The respective Atlantic and Indian inflows show a similar peak-to-peak variation of ∼2000 Gt per month, reaching their respective maxima in August and May. The Indian maximum seems to be related to a local maximum of the Pacific outflow. The annual maxima of the net WT of the four basins are reached between 3 August and 9 September, although the annual signals of the Pacific and Indian oceans are almost triple those from Arctic and Atlantic oceans (Table 1 and Fig. 6b).

Figure 7 Pacific outflow and climatic indices for ENSO, AMO, AO, and AAO. (a) The time series of the Pacific outflow is de-trended and de-seasonalized. All time series are normalized to have unit variance. Values in parenthesis are the correlation coefficients between the corresponding climatic index and the Pacific outflow. (b) Correlation coefficients between de-trended and de-seasonalized WT components of different regions and the climatic indices.

3.3 Interannual variability

Interannually, the Pacific outflow shows remarkable variability, mainly produced by P on the continents, which is inherited by R, and P−E in the oceans (Fig. 2). For example, the Pacific outflow shows a maximum of around 1372 Gt per month in 2009 that matches with a P−E maximum in the Pacific, a P−E minimum in the AIA oceans, and P minima in the continental basins draining to both the Pacific and AIA oceans. The opposite behaviour, a minimum of around 939 Gt per month, is observed in 2010. The difference, 433 Gt per month, is comparable to the discharge from the Amazon (Lorenz et al., 2014). In the tropical Pacific, the El Niño–Southern Oscillation (ENSO) is the strong recurring climate pattern involving changes in the temperature of seawater and air pressure in the tropical Pacific Ocean. ENSO had a mild El Niño phase in 2009 followed by a strong La Niña phase in 2010, which may be related to the interannual variability of the Pacific outflow. To elucidate this, we conduct a correlation study of the interannual Pacific outflow with respect to the major climate oscillations in the Earth's atmosphere and ocean: ENSO, the Atlantic Multi-decadal Oscillation (AMO), the Antarctic Oscillation (AAO), and the Arctic Oscillation (AO). The climatic oscillation is represented by monthly time series of its indices, which are non-dimensional functions of time derived from relevant meteorological observations; their values indicate the polarity and strength of the oscillation at a given epoch. The ENSO oscillations are measured here with the Southern Oscillation index (SOI), which represents the sea level pressure differences between Tahiti and Darwin, Australia. The AMO is a coherent mode of natural variability based upon the average anomalies of sea surface temperatures, with the AMO index reflecting the non-secular multi-decadal sea surface temperature pattern variability in the North Atlantic basin. The AAO describes the intensity of the westerly wind belt surrounding the Antarctic, quantified by the AAO index, which is the leading principal component of the 700 hPa atmospheric geopotential height anomalies poleward of 20^∘ S. The AO is to be interpreted as the surface signature of modulations in the strength of the polar vortex aloft in the Arctic, and the AO index is constructed by projecting the 1000 hPa height anomalies poleward of 20^∘ N. Figure 7a shows all of the indices with the amplitudes normalized to 1 standard deviation, as well as the de-trended, de-seasonalized, standard-deviation-normalized Pacific outflow. The correlation analysis between the Pacific outflow and the SOI shows no overall correlation (Pearson coefficient of 0.03) during the whole period, meaning that the influence of ENSO on the Pacific outflow may be restricted to the strong ENSO phases, such as in 2009 and 2010. A similar lack of correlation (lower than 0.1) is observed with respect to the AMO, AAO, and AO.

To explore this lack of correlation, we have estimated the correlation coefficient between each climatic index and each WT component (Fig. 7b). All of them are lower than 0.3 except for six cases in two regions. In the Arctic, P and P−E in the drainage basins of the Arctic show a correlation of ∼0.5 with the AO. This correlation is natural, as the Arctic is the area that influences the AO. The other region is the Pacific, where, as expected, the SOI shows a correlation of around 0.5 with P, P−E, and R in the drainage basins, and of around −0.4 with P in the ocean. However, this individual correlation does not extend to the Pacific outflow. In order to understand why this is the case, it is convenient to express the N component of the water transport as a function of P−E and dW. According to Eqs. (1) and (2), we have

$$\begin{array}{}\text{(3)}& \begin{array}{rl}N& =-(P-E{)}_{\mathrm{ocean}}-R+d{W}_{\mathrm{ocean}}\\ & =-\underset{X_{\mathrm{1}}}{\underbrace{(P-E{)}_{\mathrm{ocean}}}}-\underset{X_{\mathrm{2}}}{\underbrace{(P-E{)}_{\mathrm{land}}}}+\underset{X_{\mathrm{3}}}{\underbrace{d{W}_{\mathrm{land}}}}+\underset{X_{\mathrm{4}}}{\underbrace{d{W}_{\mathrm{ocean}}}}.\end{array}\end{array}$$

The correlation between N and a given index can then be express as follows:

$$\begin{array}{}\text{(4)}& \mathrm{corr}(N,\mathrm{SOI})=\sum _{i=\mathrm{1}}^{\mathrm{4}}\mathrm{corr}\left(X_i,\mathrm{SOI}\right)\cdot {\displaystyle \frac{\mathrm{SD}\left(X_i\right)}{\mathrm{SD}\left(N\right)}}\end{array}$$

where corr denotes the correlation coefficient, and SD stands for standard deviation. As shown in Eq. (4), the correlation between N and a given index is a linear combination of the correlation between each component and the index. The coefficients of the linear combination, SD(X_i)∕SD(N) are proportional to the standard deviation of each component. The components of Eq. (4) for the Pacific outflow and the SOI index are shown in Table 2. Despite the fact that some of the individual components exhibit significant correlations with SOI (in particular P−E in land and ocean), when they are combined with the corresponding coefficients, their effects cancel out, yielding a negligible correlation between water transport and SOI (below 0.03 in magnitude). Note that Table 3 also provides some insights into the causes of the interannual variability of Pacific Ocean outflow. The largest standard deviation of P−E and dW in the ocean suggests that these two components might drive the interannual variability of the Pacific Ocean outflow. This is confirmed by a correlation analysis. The correlation between N and (P−E)_ocean is −0.70. The correlation between N and dW_ocean is 0.84. The correlation between N and the corresponding land components is below 0.18. In all cases, the corresponding time series have been de-trended and de-seasonalized prior to the evaluation of the correlation.

Another possible reason for the lack of correlation resides in the definition of the study regions, for which the presence of subregions with a positive and negative influence from an index results in an overall negligible or attenuated influence of the index in the overall region. For example, a positive phase of the AMO is related to an increase in P in western Europe (Sutton and Hodson, 2005) and the Sahel (Folland et al., 1986; Knight et al., 2006; Zhang and Delworth, 2006; Ting et al., 2009), but to a decrease in P in the US (Enfield et al., 2001; Sutton and Hodson, 2005) and northeast Brazil (Knight et al., 2006; Zhang and Delworth, 2006). All of these regions are included in the Atlantic drainage basin, and the influence of a positive phase of the AMO is then attenuated.

4 Comparison with other datasets

In this section, we perform a comparison using alternative datasets. In particular, we implement the following changes:

1. the CSR GRACE mascon solution is replaced by the JPL GRACE mascon solution provided by the NASA Jet Propulsion Laboratory (Watkins et al., 2015; Wiese et al., 2019). Similarly to the CSR data, JPL data are corrected for GIA effects, the C₂₀ Stoke coefficients are replaced by a solution from SLR, and data are reduced to 1^∘ regular grids from 0.5^∘ regular grids. In addition, we applied the degree-0 Stoke coefficients correction. However, CSR and JPL mascon solutions are not directly comparable. The main reason for this is that an estimate of degree-1 coefficients has been added to JPL mascon solutions, and the GAD product has not been added back. The corrections applied by JPL are not supplied separately; thus, we cannot do or undo any of the corrections to process the JPL data as we did with CSR data. In particular, the GAD product is not available for JPL. In any case, the JPL solution is useful here, as it provides an error estimate of the mascon solution that can be propagated to obtain confidence intervals of N, which are independent of those estimated with the bootstrap analysis. Table 3 shows the CI₉₅ of the mean values of the N component estimated from error propagation and bootstrap analysis for the different ocean basins. In all cases, it is observed that the CI₉₅ values from error propagation are included in those from bootstrap analysis, meaning that the CI₉₅ values from the bootstrap analysis are a conservative estimate of the error. The JPL propagated error can be expected to be similar to that propagated from CSR error estimates (which are not available); thus, we can assume that the reported CI₉₅ values for N calculated from CSR data are a conservative estimate. In addition, comparing Tables 1 and 3, it is observed that the mean values of N are quite similar and that the CI₉₅ values largely overlap. Regarding the time variability, the values of the N component from the CSR and JPL mascon solutions show Pearson correlation coefficients greater than 0.85 (p value < 10⁻³), except for the Atlantic (0.70). Thus, despite the different processing methods used for CSR and JPL data, the reported analysis for the N component is robust with respect to the choice of GRACE datasets.

2. ERA5 P and E data are replaced by several datasets for comparison purposes. The objective is not to be exhaustive in the selection, but rather to show that the reported features of the N component are quite robust with respect to the choice of the P and E datasets. The datasets considered are:

   • i.

     continental P from GPCC (Schneider et al., 2011), GPCP (Adler et al., 2018), CMAP (Xie and Arkin, 1997), UDel (Willmott and Matsuura, 2001), and GLDAS/Noah (Rodell et al., 2004; Beaudoing and Rodell, 2016);

   • ii.

     ocean P from GPCP and CMAP;

   • iii.

     continental E from GLEAM (Miralles et al., 2011; Martens et al., 2017) and GLDAS/Noah;

   • iv.

     ocean E from OAFlux (Yu et al., 2008) and HOAPS/CM SAF (Schulz et al., 2009).

The Pacific outflow is estimated using the 162 possible combinations of P and E, including ERA5. The time period used is from 2003 to 2016, except for HOAPS/CM SAF and GPCP, which span from 2003 to December 2014 and October 2015 respectively. The degree-0 corrections in the GRACE data is made for each combination. Note that only ERA5 includes P and E for both continents and oceans. All grids have been homogenized to 1^∘ regular grids. The main concern here is the heterogeneity of the spatial coverage among datasets. To make the results comparable among datasets, the computations are restricted to the common grid points, which do not cover the entire Earth (Fig. 8a). However, in spite of the fact that the principle of water-mass conservation is not fulfilled due to the partial coverage, the Pacific outflow obtained in the common grid points from ERA5 (black line in Fig. 8b) is in quite good agreement with the same signal obtained with global coverage (the red line in Fig. 3, which is also reported as the red line in Fig. 8b). The Pearson correlation coefficient between the two signals is 0.994 (p values < 10⁻³) with an average difference of around 50 Gt per month. In general, the Pacific outflows estimated from all of the P and E dataset combinations show qualitatively the same signal as the one reported in Fig. 3. For each of the 162 estimates of the Pacific outflows corresponding to the possible P and E dataset combinations, we evaluated the average outflow (over the period of study), which is 968 Gt per month (SD of 489), and the correlation with the Pacific outflows in Fig. 3, which is 0.82 (SD of 0.06; p values < 10⁻³).

Figure 8 Monthly time series of (the opposite of) the Pacific outflow estimated from 162 combinations of P and E datasets. (a) Spatial coverage common to all datasets. (b) Pacific outflows: grey thin curves are the 162 Pacific outflows estimated in the grid points common to all datasets (no global coverage); black and red curves are based on ERA5 P and E and are obtained using either only the grid points common to all datasets (black curve) or global coverage (red curve). Note that the red curve is the same as in Fig. 3.

Table 2 Correlation coefficients between the SOI and the de-trended and de-seasonalized WT components involved in the estimate of the Pacific outflow according to Eqs. (3) and (4).

Table 3 Mean of the N component as estimated from the JPL mascon solution for different ocean basins according to Eq. (2). CI₉₅ values are estimated as the propagation of mascon errors provided by JPL and from bootstrap analysis. Units are gigatonnes per month.

These experiments show that the reported net WT values are physically consistent among datasets, at least qualitatively.

5 Discussion and conclusions

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 and 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-surficial changes would masquerade as WT variations. However, they are not critical as far as our non-secular analysis is concerned.

We use the proposed methodology to estimate the net WT and exchanges among the Pacific, Atlantic, Indian, and Arctic oceans, for the period from 2003 to 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 climatology is such that the Pacific Ocean loses water at a rate from near zero up to the peak of 2000 Gt per month during the boreal summer, which coincides with the maximum of the global atmosphere water content. On top of the climatology, the interannual Pacific annual mean water loss varies significantly between ∼950 and ∼1450 Gt per month during the study period, but this is seemingly uncorrelated with ENSO.

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 see-saw WT between the Pacific and the Atlantic and Indian oceans has been reported (Chambers and Willis, 2009; Wouters et al., 2014). In those studies, the $P-E+R$ term in Eq. (2) in both the Pacific and the Atlantic and Indian oceans was approximated by that from the global ocean mean. However, the mean freshwater flux in the Pacific (1261 Gt per month) does not really match with that in the Atlantic and Indian oceans (−1837 Gt per month), meaning that the approximation was quite poor; hence, the N term was not properly estimated in these studies (see the Appendix for further discussion).

Differences in freshwater fluxes between the Pacific and Atlantic oceans produce salinity contrasts and, in turn, differences in 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 and Cessi, 2017), the meridional span of the African and American continents (Nilsson et al., 2013; Cessi and Jones, 2017), and the salty WT from the Indian Ocean to the Atlantic (Gordon, 1986; Marsh et al., 2007). The AMOC is also influenced by WT through the Bering Strait (Reason and Power, 1994; Goosse et al., 1997; Wadley and Bigg, 2002) as well as by the surface processes of temperature, precipitation, and evaporation at low latitude in the Pacific and Indian oceans (Newsom and Thompson, 2018). The relative importance of the multiple drivers influencing the AMOC is an open problem (Ferreira et al., 2018). The net WT estimated here provides information on differences between oceanic inflows and outflows, which can be useful to help clarify this problem.

Net WT in the open oceans can alternatively be estimated using global ocean models, which simulate observational data based on physical principles. However, these models are not necessarily sensitive to the WT, especially given the data types and the geographical and topographical resolutions involved in the models. Knowledge about three-dimensional global ocean circulation could also elucidate the net WT. However, the small ratio between the net and the total WT hinders the estimation of the former from the latter.

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 A: Apparent net mass exchange between the Pacific and the Atlantic and Indian oceans

We shall show here that the net water-mass exchange between the Pacific and the Atlantic and Indian oceans reported by Chambers and Willis (2009) was a mathematical artefact. Their Eq. (2) approximated the freshwater flux, i.e. $P-E+R$, of the Pacific (Pcf) and the Atlantic and Indian (AI) oceans using the global ocean (GO) mean. However, from Figs. 2 and 4, we get very different results, ($P-E+R{)}_{\mathrm{Pcf}}=\mathrm{1261}$ and ($P-E+R{)}_{\mathrm{AI}}=-\mathrm{1837}$ Gt per month, meaning that the approximation in Chambers and Willis (2009) and, hence, their resultant estimates of WT are rather poor. In addition, using their approximation, an apparent net mass exchange will always arise, as

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}(P-E+R{)}_{\mathrm{GO}}& =\sum _{x\in \mathrm{GO}}{\displaystyle \frac{(P-E+R{)}_x\cdot \mathrm{Area}(x)}{\mathrm{Area}\left(\mathrm{GO}\right)}}\\ & =\sum _{x\in \mathrm{Pcf}}{\displaystyle \frac{(P-E+R{)}_x\cdot \mathrm{Area}(x)}{\mathrm{Area}\left(\mathrm{GO}\right)}}\\ & +\sum _{x\in \mathrm{AI}}{\displaystyle \frac{(P-E+R{)}_x\cdot \mathrm{Area}(x)}{\mathrm{Area}\left(\mathrm{GO}\right)}}\\ & =\sum _{x\in \mathrm{Pcf}}{\displaystyle \frac{(P-E+R{)}_x\cdot \mathrm{Area}(x)}{\mathrm{Area}\left(\mathrm{Pcf}\right)}}\\ & \cdot {\displaystyle \frac{\mathrm{Area}\left(\mathrm{Pcf}\right)}{\mathrm{Area}\left(\mathrm{GO}\right)}}+\sum _{x\in \mathrm{Atl}/\mathrm{Indian}}\\ & {\displaystyle \frac{(P-E+R{)}_x\cdot \mathrm{Area}(x)}{\mathrm{Area}\left(\mathrm{AI}\right)}}\cdot {\displaystyle \frac{\mathrm{Area}\left(\mathrm{AI}\right)}{\mathrm{Area}\left(\mathrm{GO}\right)}}\\ & \approx (P-E+R{)}_{\mathrm{Pcf}}\cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\\ & +(P-E+R{)}_{\mathrm{AI}}\cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},\end{array}\end{array}$$

where x are disjoint grid cells in the ocean basins; the areas of the Pcf, AI, and GO are around 177, 160, and 351×10⁶ km²; and the ratios 177∕351 and 160∕351 have been approximated by 1∕2. Then, multiplying by 2 and rearranging the equation, we get

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}(P-E& +R{)}_{\mathrm{Pcf}}-(P-E+R{)}_{\mathrm{GO}}\\ & \approx -\left[(P-E+R{)}_{\mathrm{AI}}-(P-E+R{)}_{\mathrm{GO}}\right].\end{array}\end{array}$$

Thus, wherever the signal is in the Pacific and the Atlantic and 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 an exchange does not exist.