2.1 Evaporation data records

The liquid-water-equivalent evaporation rate, E, is calculated from the latent heat flux Q_l as follows:

where L_E is latent heat of evaporation of water. The latent heat flux, in turn, is parameterized according to the bulk flux algorithm (based on the Monin–Obukhov similarity theory representation of fluxes in terms of mean quantities):

with being ρ the density of air; C_E the coefficient of turbulent exchange; u the wind speed at 10 m height relative to the ocean surface current speed; and q_s and q_a the specific humidity at the sea surface and at 10 m height, respectively. Whereas q_a and u are derived from satellite observations of BT, ρ, q_s, and L_E are derived from their dependences on SST and/or air temperature. The turbulent exchange coefficient C_E is obtained from the Coupled Ocean-Atmosphere Response Experiment (COARE) version 3.0 algorithm (Fairall et al., 1996, 2003). The algorithm iteratively estimates stability-dependent scaling parameters and wind gustiness to account for sub-scale variability.

Most of the data sets used here do not explicitly contain E; therefore, we calculated those from monthly means of Q_l and SST using Eq. (5) and L_E (in J kg⁻¹) given by (Henderson-Sellers, 1984)

where T_s is SST in kelvin. The slight difference with the definition of L_E used in the COARE-3.0 algorithm causes negligible differences of 0.03–0.04 % for T_s between 278 and 298 K.

The BT observations common to satellite-based retrievals of ocean turbulent fluxes come from the Special Microwave Imager (SSM/I; Hollinger et al., 1990) and Special Microwave Imager/Sounder (SSMIS; Liman et al., 2008) instruments on the Defense Meteorological Satellite Program (DMSP) platforms F08–F18. These data were corrected and intercalibrated using various approaches to create FCDRs, stable fundamental climate data records (see, e.g., Wentz et al., 2013; Sapiano et al., 2013; Berg et al., 2018; Fennig et al., 2020), which then serve as input to various satellite retrievals. Slight differences in calibration approaches lead to differences in FCDRs that propagate into the retrieved data. Issues with sensor stability, especially with SSM/I and SSMIS sensors, usually express themselves as slow drifts or sudden jumps of the global mean.

2.2 Precipitation data records

Microwave-based retrievals of precipitation are based on the interaction of liquid or solid hydrometeors with the upwelling radiation field. In HOAPS-4.0, P is determined by a neural network retrieval trained on profiles from an ERA-Interim climatology (Andersson et al., 2010). The training data set consists of 1 month (August 2004) of assimilated SSM/I BTs and the corresponding ERA-Interim P (Bauer et al., 2006).

There are a multitude of global precipitation products in existence (see, e.g.,  Kidd and Huffman, 2011; Tapiador et al., 2017), but for this study we selected GPCP as the P data set with which to calculate E−P (except for the HOAPS product, which makes use of its own P data) because it is generally regarded as the data set that performs best globally. Moreover, J-OFURO also makes use of GPCP P to determine E−P (Tomita et al., 2019).

The Global Precipitation Climatology Project–One-Degree Daily data set (GPCP-1DD; denoted GPCP hereafter) contains P estimated from a combination of data from ground-based rain gauges and satellites – the latter including near-infrared, passive, and active microwave observations (Huffman et al., 2001). Daily global precipitation rates are provided by GPCP-1DD at 1^∘ resolution for the time range 1996–2017. We calculate monthly mean P from version 1.3 GPCP-1DD (GPCP, 2018), because the spatial resolution of the monthly product is not sufficient for our purposes. These data were obtained from https://rda.ucar.edu/datasets/ds728.5/ (last access: January 2021).

2.3 Errors, biases, and uncertainty

Four out of seven data sets analyzed here contain explicit information on uncertainty. HOAPS contains estimates of random and systematic bias errors (Kinzel et al., 2016; Liepert and Previdi, 2012). The errors in E were obtained by separating biases of HOAPS Level 2 E with respect to collocated in situ ship-based data into equally populated E, u, T_s, and W bins. The mean and standard deviation of the biases are assumed to represent the systematic and random components of the 2σ uncertainty range, respectively, which is probably a conservative estimate. When the approach is taken of determining uncertainty ranges as a function of turbulent flux parameters, these can also be assigned to times and regions not covered by the ship-based reference data set (Liepert and Previdi, 2012). For the current study, we calculated the mean uncertainty by averaging the systematic uncertainty component. The random component is negligible when averaging long time series. The HOAPS P data set does not contain uncertainty information; instead, a constant relative 1σ uncertainty range of 13 % was assumed, based on a comparison with ship-based in situ data (Burdanowitz, 2017). The total E−P uncertainty was determined by error propagation.

Bias errors given in the OAFlux data set were computed based on the uncertainty ranges of individual input data sets, assuming no correlation between uncertainties from different data sets (Yu et al., 2008). Like for HOAPS uncertainty ranges, the OAFlux bias error was simply averaged for our investigations.

Uncertainties in SEAFLUX arise both from comparisons of the individual retrieval errors (e.g., wind speed, humidity, air temperature) evaluated against quality-controlled buoy archives and from errors arising as a result of gap-filling through application of a Kalman smoother. Individual retrievals were generally found to be unbiased globally, but some conditional biases likely remain. The total uncertainty is a measure of the reduction in retrieval uncertainties through combination of multiple sensors at each location and time and increases in uncertainty related to sampling inhomogeneities. As the length of time grows between any given time and the previous or next observation, the sampling uncertainty increases. Thus the SEAFLUX uncertainties generally capture random retrieval uncertainties and sampling uncertainty but do not contain conditional systematic errors as developed for HOAPS. However, we note that the retrieval error itself does likely contain some components of conditional systematic biases even though the unconditional biases remain small.

In contrast to the monthly GPCP product, GPCP-1DD Version 1.3 does not provide explicit uncertainty estimates; hence here we assume a constant relative 1σ uncertainty range of 8 %. This is the estimated bias error for GPCP data over the tropical oceans (Adler et al., 2012), which is where most of the P signal originates. Over the global oceans, the bias error was estimated at 10 %, but Adler et al. (2012) considered this an upper bound.

In contrast to uncertainty ranges estimated by comparing with other (e.g., in situ) data sets, the uncertainty of ERA5 data is described by the standard deviation and the range of the ensemble, consisting of 10 separate reanalysis runs (Hersbach et al., 2020). We determined these statistics after averaging of the data: first, the mean (e.g., global monthly mean) of each individual ensemble member was calculated, and then standard deviation and range were determined. Note that ERA5 ensemble statistics should be interpreted in a relative sense (i.e., ensemble spread is larger where uncertainty is higher), as the numerical values are overconfident (ECMWF, 2020).

4.1 Freshwater flux climatology

Freshwater flux climatologies obtained from 17 years of data (1997–2013) were determined from satellite ensemble median (SEM) and ERA5 data. They are shown in Fig. 1a and b, to illustrate the overall spatial distribution of mean E−P. The chosen time range is the largest common time range of the data sets used in this study. Note that ERA5 data were matched to satellite data coverage.

Figure 1Satellite ensemble median (SEM) and ERA5 climatologies (1997–2013) of freshwater flux (a, b) and evaporation (c, d), and GPCP and ERA5 precipitation (e, f). ERA5 data coverage was reduced to match satellite data, and data over land were discarded from panels (e) and (f). See the text for details.

Regions where mean P>E are dominated by atmospheric freshwater outflux (into the ocean), appear in blue in Fig. 1a and b, and are concentrated at the Intertropical Convergence Zone (ITCZ) and the Pacific warm pool. In the subtropics, E generally outweighs P. At higher latitudes P and E are approximately equal, but with a tendency to $E-P<\mathrm{0}$. Comparison of panels c–f with a and b shows that the E−P pattern is mainly determined by P in the tropical and high-latitude regions but determined by E in the subtropical regions. The agreement between SEM and ERA5 E−P climatologies is good, yet some systematic differences can be observed. Due to higher P in the ITCZ, ERA5 shows more negative E−P there. Conversely, the overall higher E level in ERA5 causes E−P values larger than those found for SEM over most of the global oceans. Excessive E was also found to produce high E−P in ERA-Interim (Brown and Kummerow, 2014).

Figure 2Difference maps of HOAPS (a), OAFlux-G (b), and SEAFLUX-G (c) climatological mean E−P minus the corresponding collocated ERA5 climatology (1997–2013). HOAPS (d), OAFlux-G (e), and SEAFLUX-G (f) climatological mean 1σ uncertainty. White lines in the left panels enclose regions where the difference with ERA5 E−P exceeds the 2σ uncertainty range.

The deviations are more apparent when climatological differences are analyzed. For this comparison we select ERA5 as a reference due to its spatiotemporal completeness and because it is the only “other” data set (i.e., not satellite data), keeping in mind that ERA5 data very likely also have inaccuracies and/or biases. Figure 2 shows climatological difference plots of HOAPS (upper panel), OAFlux-G (middle panel), and SEAFLUX-G (lower panel) with collocated ERA5 data. Although HOAPS differences with ERA5 appear larger to the eye, root mean squared differences are 0.6 mm d⁻¹ for each of the three comparisons: 0.60 mm d⁻¹ for HOAPS, 0.58 mm d⁻¹ for SEAFLUX-G, and 0.57 mm d⁻¹ for OAFlux-G. As already seen in Fig. 1, differences are not homogeneously distributed over the globe. The HOAPS difference plot is characterized by an alternating pattern of positive and negative deviations. Stronger HOAPS E in the subtropical central north and eastern South Pacific produces elevated E−P compared to ERA5. In contrast, elevated ERA5 E over the East China Sea combines with smaller ERA5 P in the region, resulting in higher ERA5 E−P. The positive bands on either side of the Equator are due to higher HOAPS E, whereas the negative E−P differences at the Equator are due to smaller HOAPS P. The negative deviations to the east and west of Australia are also due to differences in P, whereas the deviations at latitudes $>\mathrm{40}{}^{\circ }$ S are due in equal parts to E and P. The differences between OAFlux-G and ERA5 are mainly due to P, apart from the regions in the subtropical Pacific and Atlantic Oceans, where OAFlux E is smaller than ERA5 E. SEAFLUX-G shows slightly larger differences with ERA5. In the band within 30^∘ of the Equator, SEAFLUX yields higher E than ERA5 (and OAFlux) in most of the Pacific and Atlantic Ocean, except in the upwelling regions on the west coasts of Africa and the Americas. The difference plots of J-OFURO and IFREMER-G with ERA5 are not shown here but are very similar to the lower left panel because the differences in P between GPCP and ERA5 are larger than differences in E in most regions. All plots, including difference climatologies of E and P, can be found in the Appendix, Fig. A1.

To investigate where the differences are significant, the right column of Fig. 2 presents the 1σ uncertainty range from HOAPS (upper panel), OAFlux-G (middle panel), and SEAFLUX-G (lower panel). Moreover, regions where the difference between satellite E−P and ERA5 E−P is greater than the 2σ uncertainty range are enclosed by white contour lines in the left panels. The ERA5 E−P uncertainty shows a pattern similar to that of OAFlux-G but is a factor of 10 smaller than the uncertainties estimated for satellite data and therefore adds a negligible component to the total uncertainty estimate. The HOAPS uncertainty range is larger than HOAPS-ERA5 E−P differences over most of the globe. This is mainly due to P, for which we assumed 13 % uncertainty. The deviations >1 mm d⁻¹ in the oceans' desert regions (off the west coasts of Peru and southern Africa) and in the higher latitudes are clearly outside the 2σ uncertainty ranges. In contrast, OAFlux-G E−P deviations are larger than the estimated 2σ uncertainties in the ITCZ, on the west coasts of the Pacific and Atlantic Ocean, in the Arabian Sea, and at southern high latitudes. Again, the uncertainty range is mainly given by P, for which we assumed a relative uncertainty of 8 %. Due to the small uncertainty estimates in SEAFLUX, all of the larger differences with ERA5 in the Atlantic and Pacific Ocean are significant.

4.2 Intercomparison of freshwater flux over ocean: global means

Monthly mean E, P, and E−P of six (or three) data sets were collocated (see Sect. 3) and averaged over the global oceans (80^∘ S–80^∘ N). Climatological seasonal cycles were determined for the overlapping time range (1997–2013) and are shown in Fig. 3a–c. HOAPS, ERA5, OAFlux, SEAFLUX, and GPCP uncertainty ranges are presented in the boxes attached to the right of panels a–c. Dots show the climatological mean value, and error bars indicate the associated 1σ uncertainty. Subtracting the seasonal cycle from the respective monthly mean time series yields global ocean anomalies of E−P, E, and P, which are presented as 3-month running means in panels d–f. Seasonal and interannual variability are on the same order of magnitude, which can be seen by comparing the left panels with those on the right (the y axis spans 1 mm d⁻¹ in all panels).

Figure 3Climatological (1997–2013) seasonal cycle of global ocean mean evaporation rate (a), precipitation rate (b), and freshwater flux (c). HOAPS, ERA5, OAFlux, SEAFLUX, and GPCP mean values and associated 1σ uncertainty ranges are shown in the boxes to the right of the panels. Monthly mean anomaly (with respect to the climatological seasonal cycle depicted on the left) over the global oceans (80^∘ S–80^∘ N) of evaporation rate (d), precipitation rate (e), and freshwater flux (f). The anomaly data are smoothed using a 3-month running mean. Panel (e) additionally displays the Niño3.4 index shifted by +3 months (right y axis). The legend shows the correlation coefficient of the Niño3.4 index with P anomalies and the time lag of highest correlation (Δt in months). Ticks on the time axis mark January of the indicated year.

Download

There are substantial deviations between E, P, and E−P data. Fig. 3a shows that a difference of about 0.2 mm d⁻¹ is found between OAFlux and J-OFURO E. An additional discrepancy of 0.2 mm d⁻¹ exists between J-OFURO and ERA5. E data from HOAPS, IFREMER, and OAFlux are much closer to each other: satellite-based E all falls within the OAFlux uncertainty range (red error bars), whereas the ERA5 climatological mean E does not fall within the larger HOAPS uncertainty range. The HOAPS uncertainty range is much larger than the seasonal variation, which indicates that it is likely overestimated, which may be due to the assumption of 100 % covariance for systematic uncertainty.

Fig. 3b shows that the seasonal cycle of global ocean mean P is shallow, and the two satellite-based data sets agree within the GPCP uncertainty for 10 months of the year. Like for E, we find substantial differences among the three P data sets: there is a deviation of about −0.1 mm d⁻¹ between HOAPS and GPCP, and ERA5 shows values that are about 0.25 mm d⁻¹ higher than GPCP, which was also found by Hersbach et al. (2020). These differences can, in part, be explained by differences in P frequency distributions and, in particular, by the fraction of rain occurrences, which is much lower in HOAPS than in GPCP or ERA5. This will be discussed in Sect. 5. Since in this paper the focus is on the intercomparison of E−P (not specific E or P algorithm issues), we only describe the observed differences between P (and E) data sets to obtain a better understanding of differences between E−P data.

Apart from HOAPS E−P in March–April, all satellite data sets agree on phase and amplitude of the E−P seasonal cycle (Fig. 3c). ERA5 shows hardly any dependence on season, as the magnitude of the summer maximum is smaller in ERA5 due to the relatively larger summer P maximum. The monthly and interannual variability of ERA5 E−P is, like the seasonal cycle, of smaller amplitude than that of satellite data, which is caused by the high degree of coherence between E and P, and will be discussed in more detail in Sect. 4.5. Because, compared to satellite data, ERA5 E and P are biased high by about the same amount, E−P is close to the satellite data. HOAPS yields the highest E−P due to its low mean P. All E−P data are contained within the HOAPS and OAFlux uncertainty ranges.

The E anomalies in Fig. 3d display a high degree of correlation on a monthly timescale. On a multi-annual scale all data sets show some degree of variability, which is most likely linked to sensor and intercalibration issues (e.g., Robertson et al., 2020), and the variability is not consistent. For example, the slow, decadal-scale oscillation observed in HOAPS and IFREMER appears to be in anti-phase compared to OAFlux. The three P data sets yield interannual variations with amplitudes that are similar in amplitude to those found for E and show a high degree of correspondence in their monthly and interannual variability – apart from the stronger dependence of HOAPS on ENSO (El Niño–Southern Oscillation) phase. This is a known characteristic of HOAPS data (see, e.g., Andersson et al., 2011; Masunaga et al., 2019) and is most apparent in panel e, where the Niño 3.4 SST index (Trenberth and Stepaniak, 2001) is plotted in gray bars along with P anomalies: HOAPS P correlates with Niño 3.4 if a lag of 3 months is taken into account (R²=0.73). Apparent agreement is found among all E−P anomalies (panel f) – again apart from the ENSO-related deviations found in HOAPS P. The agreement among E−P anomalies is best in the “quiet” ENSO years (2001–2005), but this is probably a coincidence as the spread in E−P in other years is mainly due to differences in E and not in P. Note that differences between J-OFURO-G, IFREMER-G, SEAFLUX-G, and OAFlux-G are due to differences in E, as in all cases GPCP P was used for the calculation of E−P.

4.3 Intercomparison of freshwater flux over ocean: time series on regional scales

In this section, we investigate the temporal correlation of water cycle components on regional scales. This approach will help to elucidate differences between the various data sets by uncovering in which regions the differences are particularly large (or small). As a reference for the E and E−P comparisons, we use SEM, a data set determined by the statistical median of all satellite data sets. Since we use only two satellite P data sets, GPCP is selected as a reference for the P comparison. We determine correlation coefficient, slope, and intercept of the linear regression ($y=ax+b$) between $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ monthly means (not anomalies) of each data set, y, and the reference, x, to examine where estimates are most consistent.

Figure 4Correlation, slope, and intercept of the linear regression of monthly mean E from (top to bottom) HOAPS, J-OFURO, IFREMER, OAFlux, SEAFLUX, and ERA5 with satellite ensemble median (SEM) monthly mean E (1997–2013).

The results are shown for all six E data sets in Fig. 4, where the left column displays the correlation coefficient. In the top row, HOAPS yields R²>0.75 over most of the globe, with some notable exceptions in the ITCZ and at the Peruvian coast. The other satellite data yield higher correlation coefficients. The correlation pattern of ERA5 with SEM is similar to that found for HOAPS, although the tropical areas with R²<0.75 are not at the same locations. The highest overall correlation with SEM is found for J-OFURO and SEAFLUX, with R² exceeding 0.75 essentially everywhere.

The middle panels of Fig. 4 display the slope of the linear regression. A slope greater (smaller) than 1 implies an overestimation (underestimation), particularly of large values, compared to SEM. HOAPS overestimates E in the tropics, except in an area in the eastern Pacific at 0–5^∘ N, where a<1. J-OFURO, IFREMER, and OAFlux each yield slopes <1 within 30^∘ of the Equator and slopes close to 1 everywhere else (apart from the band with a<1 seen in IFREMER at high southern latitudes). Of those three E data sets, OAFlux displays the largest deviations from a=1. In contrast, SEAFLUX yields slopes close to unity over the whole globe. An inhomogeneous pattern is found for ERA5, but the slope is generally close to 1. A small region in the tropical Atlantic stands out due to its large slope, and since this is not seen in any of the satellite data sets, it must be a feature in ERA5 data.

The patterns in the middle panels are nearly all mirrored in the right panels; i.e., wherever large values are overestimated (a>1), small values are underestimated (b<0), and vice versa. All data sets thus appear to agree on intermediate values. Overall, the correspondence between E data sets is best in the subtropics, while the largest deviations appear mainly in the tropics. This is due to the frequent occurrence of weather conditions in which the moisture stratification departs substantially from typical conditions to which the retrieval algorithms of near-surface moisture are tuned. Accounting for this dependence on moisture stratification, as in the SEAFLUX and J-OFURO algorithms, improves retrieval results appreciably compared to in situ measurements (Roberts et al., 2019).

Figure 5Correlation, slope, and intercept of the linear regression of monthly mean P from HOAPS (a–c) and ERA5 (d–f) with GPCP monthly mean P (1997–2013).

Figure 5 shows the same analysis for P from HOAPS (upper panels) and ERA5 (lower panels). The correlation coefficient between HOAPS and GPCP P is >0.75 in the ITCZ and about 0.5 for most of the global oceans. In the oceans' deserts R²<0.25 is found, which is mostly due to the small dynamic range of mean P. Compared to GPCP, HOAPS underestimates P in this region, as a<1. At latitudes poleward of 50^∘ similarly small R² values are found that are due in part to the small dynamic range and in part to difficulties pertaining to the detection of snow by passive microwave instruments (Tapiador et al., 2017; Kidd and Huffman, 2011). HOAPS underestimates high P here and overestimates small P (b>0 mm d⁻¹) compared to GPCP. Very similar patterns are seen for ERA5, although in general the correlation coefficient is higher than for HOAPS. ERA5 is biased high almost everywhere compared to GPCP. Both HOAPS and ERA5 show a smaller range of P in the Southern Oceans, as the slope is less than 0.5, but the large intercept indicates an overestimation of small P compared to GPCP. The narrow band of R²<0.75 and b>1 mm d⁻¹ at the Equator is also found in both HOAPS and ERA5.

Figure 6Correlation, slope, and intercept of the linear regression of monthly mean E−P from (top to bottom) HOAPS, J-OFURO-G, IFREMER-G, OAFlux-G, SEAFLUX-G, and ERA5 with satellite ensemble median (SEM) monthly mean E−P (1997–2013).

For E−P, the results of the regression analysis are shown in Fig. 6. The highest correlation coefficients (and slopes and intercepts closest to 1 and 0 mm d⁻¹) are found among the data sets calculated with GPCP P. This shows that most of the variability in E−P is due to differences in P. Since GPCP P is used in four out of five data sets included in the SEM, those data sets show high correlations, whereas HOAPS and ERA5 yield patterns very similar to those found for the P comparison in Fig. 5. Nevertheless, both for ERA5 and HOAPS the correlation in most of the tropics is higher for E−P than for P. In summary, the correlation patterns for HOAPS and ERA5 indicate agreement on the seasonal cycle in the tropics, a result found previously by Brown and Kummerow (2014), although we find that its amplitude is reduced in the GPCP-based E−P data (Fig. 3). Less agreement is found in the Southern Oceans, where GPCP-based E−P is underestimated relative to SEM. At the midlatitudes, the regression with SEM yields slopes near 1 and intercepts close to 0 mm d⁻¹ for ERA5 and HOAPS, but the correlation is less than in the tropics, probably due to the smaller dynamic range of E−P.

In the present study, we compare satellite-based E−P with ERA5 E−P because we are also examining the separate contributions from E and P. It can, however, be argued that VIMD from reanalysis is a more reliable quantity than reanalysis E−P, since VIMD is calculated from the state variables wind and water vapor, whereas E and P are derived from model physics (e.g., Trenberth et al., 2011). We verified that in ERA5 the agreement between E−P and ∇⋅(vq) is generally good, as shown in Appendix B. Hence, changes to the plots in Fig. 6 are minor when ERA5 ∇⋅(vq) is used to calculate the regression with SEM E−P instead of ERA5 E−P, as shown in Fig. A2.

4.4 Examination of the water budget in ERA5

One way of investigating the consistency of different water cycle components is determining if the global water budget (Eq. 1) is closed. However, satellite E−P data sets are available over ocean only, so we revert to a comparison with gap-free reanalysis data. There is no internal constraint for budget closure in ERA reanalyses (Berrisford et al., 2011; Hersbach et al., 2020), and as the budget was not closed in ERA-Interim, it is worthwhile to investigate ERA5's behavior in this regard. Monthly mean total ERA5 E−P over the globe, the ocean, and land is shown in Fig. 7 in black, blue, and green, respectively. The mean values over the globe and land were scaled by their surface area relative to the ocean surface area (i.e., they were multiplied by 510∕350 and 160∕350, respectively) to obtain consistency with the over-ocean means shown in Fig. 3. The error bars on ERA5 data depict the standard deviation of the 10-member ensemble. Nearly all of the uncertainty in global mean E−P is due to the uncertainty over ocean; the error bars on the over-land E−P are smaller than the graph's line width. This is due in equal parts to E and P, which have similar ensemble standard deviations (not shown). For the time range shown in Fig. 7 global E−P is seen to oscillate around 0 mm d⁻¹, meaning that the ERA5 water budget is closed on a yearly timescale (in agreement with the findings by Hersbach et al., 2020). The seasonal cycle is mainly driven by increased evapotranspiration of vegetation on land and peaks in northern hemispheric summer due to the larger fraction of land in the Northern Hemisphere. Precipitation shows a similar seasonal cycle over land but does not completely cancel out in E−P due to a slight phase shift with respect to the E seasonal cycle (not shown).

Figure 7ERA5 monthly mean E−P over the whole globe (black), land only (green), and ocean (blue); global mean ΔW (light blue), and mean ∇⋅(vq) over land (pink) and ocean (purple). The mean values over the globe and land were scaled by their surface area relative to the ocean surface area (i.e., they were multiplied by 510∕350 and 160∕350, respectively) to obtain consistency with the over-ocean means shown in Fig. 3. Error bars represent the standard deviation within the 10-member ensemble, which is smaller than the graph's line width for E−P over land, ΔW, and ∇⋅(vq).

Download

Figure 7 shows that monthly means of global E−P and ΔW (light blue line) display a high degree of coherence, as expected from Eq. (2). This is an indication that the (atmospheric) water cycle is well represented in ERA5.

Globally, VIMD is zero, as no water vapor is transported out of (or into) the Earth system. However, we find ERA5 global total VIMD to be −0.04 mm d⁻¹: a small value within the standard deviation of the ensemble of single grid boxes, but significant and on the order of the amplitude of the seasonal cycle of net E−P on a global scale. The deviation from zero is due to the fact that VIMD is calculated in grid point space (and not in the model's spectral space), where the mathematical constraint of net zero divergence is not enforced (Paul Berrisford, personal communication, October 2020). Interestingly, VIMD over land (pink) agrees well with over-land E−P, whereas VIMD over ocean (purple line) is smaller than over-ocean E−P also by −0.04 mm d⁻¹. Based on the results of the regression analysis shown in the upper panels of Fig. A2 we speculate that discrepancies between E−P and ∇⋅(vq) over the ocean's desert regions also play a role in causing ∇⋅(vq) to be smaller over ocean than over land.

4.5 Examination of the water budget in satellite data sets

Globally, E−P is equal to ΔW (Eq. 2 and Fig. 7), and because ΔW is 2 orders of magnitude smaller than E and P, global mean E should necessarily be almost equal to global mean P. This seemingly trivial finding provides us with a tool for investigating the consistency of E and P data sets: by determining how well they correlate. For ERA5, global mean E and P yield correlation coefficients R²=0.82 and R²=0.84 for monthly and yearly means, respectively. This procedure cannot be applied to the satellite E−P data considered here, as they contain values over ocean only. Since there is a substantial seasonality in water vapor transport (Fig. 7), the correlation between ocean-only E and P is expected to be much lower. A regression of ERA5 E_ocean and P_ocean monthly means (where the subscript _ocean indicates averaging over ocean only) indeed yields only R²=0.42, or R²=0.57 for yearly means. To account for the net transport of water from ocean to land, we include ∇⋅(vq)_ocean into the analysis and, applying Eq. (4), correlate $E_{\text{ocean}}-\mathrm{\nabla }\cdot (vq{)}_{\text{ocean}}$ with P_ocean. For ERA5, the resulting R² are 0.86 for yearly and monthly means: very similar to the coefficients found for the correlation of global E with P.

We calculated correlation coefficients of the various E and P data sets used in this study, combined with ERA5 ∇⋅(vq)_ocean, and listed them in Table 3. The analyses were performed separately on (i) monthly means, primarily indicating agreement on the seasonal cycle; (ii) yearly means, a measure of consistency of interannual variability, including trends; and (iii) monthly anomalies, focused on short-term variability.

Table 3Pearson's correlation coefficient squared (R²) for monthly (mean or anomaly) or yearly global ocean mean $E_{\text{ocean}}-\mathrm{\nabla }\cdot \left(vq\right)$ vs. P_ocean, with ∇⋅(vq) data from ERA5. R² was calculated from data sets that were collocated prior to the calculation of global means.

^* Non-significant correlation coefficients (p value > 0.05).

Download Print Version | Download XLSX

The small correlation coefficients found for monthly mean satellite data in part reflect the differences in the seasonal cycles of E and P (see panels a and b of Fig. 3). But the results from the analysis of monthly anomalies (where the mean seasonal cycle was subtracted) are very similar to those found for monthly means, indicating that, compared to monthly and interannual variability, the seasonal cycle is of lesser importance on a global scale. Including the contribution of ∇⋅(vq)_ocean improves the correlation appreciably for ERA5, as mentioned above. For satellite data the correlation also improves, particularly for yearly means and monthly anomalies of IFREMER-G and J-OFURO-G (not shown). On a yearly timescale, we do not expect a high degree of correlation, as interannual variability is small and no clear trends are observed in panels d and e of Fig. 3. For ERA5, R²=0.86, but this is primarily caused by small E and P in 1997 and 1999, which is also the case for IFREMER. The correlation found for J-OFURO, R²=0.31, is the highest found among satellite data. The remaining satellite data sets are not significantly correlated on a yearly timescale (p value > 0.05). Clearly, time series longer than the 17 years investigated here would benefit the analysis of yearly mean data.

Table 4Estimates of ocean total E and P, land and ocean total E−P, net transport of water vapor, and continental runoff given in 10³ km³ yr⁻¹. The upper three rows contain results from this study, and the lower five those from earlier investigations. ERA5 estimates are calculated from ensemble mean data; the standard deviation (SD) is derived from ensemble statistics. The satellite-based data sets used in our study were averaged to obtain the mean and SD of observed (Obs.) E_ocean and P_ocean, and the range is given in the third row. Net water vapor flux divergence over land (∇⋅(vq)_land) and ocean (∇⋅(vq)_ocean) and continental runoff R are given in the last three columns. The estimates from the study by Rodell et al. (2015) are separated into observations (obs.) and model-optimized observations (opt.); see the text for details.

Download Print Version | Download XLSX

Overall, this analysis shows that satellite-based estimates of E are less consistent with satellite-based P data than ERA5 E and P. To a certain degree, this is expected, as the three variables used in the analysis come from different sources (e.g., E from IFREMER, P from GPCP, and ∇⋅(vq)_ocean from ERA5), each with its own sampling and uncertainty characteristics. Nevertheless, from the global water cycle perspective, some degree of correspondence between E_ocean and P_ocean is expected.

4.6 Estimates of global total water cycle components

From the data shown in the previous sections we calculated climatological (1997–2013) global ocean total E and P, total E−P (separated into land and ocean contributions), runoff, and net transport. The latter is equal to the total over-ocean (or over-land) ∇⋅(vq). Our results are given in the three upper rows of Table 4. For brevity, we again denote total E, P, and E−P over ocean as E_ocean, P_ocean, and (E−P)_ocean, respectively, and, similarly, over-land E−P as (E−P)_land. For better comparison with earlier estimates, values are given in units of 10³ km³ yr⁻¹.

The rows labeled “obs.” display the mean, standard deviation (SD), and range of P_ocean (from GPCP and HOAPS) and E_ocean (from HOAPS, J-OFURO, IFREMER, SEAFLUX, and OAFlux). Various estimates of global totals of water cycle components can be found in the literature – for example in Oki and Kanae (2006), Trenberth and Asrar (2014), Rodell et al. (2015), and Allan et al. (2020) – and are shown in the lower half of Table 4.

The largest spread among water cycle components is found for E and P over ocean, both in an absolute and a relative sense, each with a range spanning about 100×10³ km³ yr⁻¹, or 20–25 %. The relative spread in R is similar but is a factor of 10 smaller in absolute values. In their study, Rodell et al. (2015) estimated E_ocean from observations at 410×10³ km³ yr⁻¹ (corresponding to 3.21 mm d⁻¹), in agreement with Trenberth and Asrar (2014), but found a value of 450×10³ km³ yr⁻¹ by applying an algorithm that optimized all water cycle components to achieve water and energy budget closure. The algorithm caused a concurrent increase in P_ocean from the observed 385×10³ to 403×10³ km³ yr⁻¹. The global total fluxes estimated by Allan et al. (2020) derive from Rodell et al. (2015), but following the recommendation by Stephens et al. (2012), E_ocean and P_ocean were both increased by 30×10³ km³ yr⁻¹ to improve the agreement with energy constraints, yet keeping land–ocean fluxes constant. These increases are larger than the $\pm \mathrm{22}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹ uncertainty on E_ocean and P_ocean estimated by Rodell et al. (2015) based on the optimized method, and so a more modest increase of about 20×10³ km³ yr⁻¹ may be more appropriate. This would produce fluxes of $E_{\text{ocean}}=\mathrm{470}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹ and $P_{\text{ocean}}=\mathrm{424}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹, which are quite close to ERA5 estimates (Richard Allan, personal communication, October 2020).

In our study, we find a large range of E_ocean: HOAPS yields $\mathrm{397}\pm \mathrm{96}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹, OAFlux $\mathrm{414}\pm \mathrm{37}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹, IFREMER 418×10³ km³ yr⁻¹, SEAFLUX $\mathrm{444}\pm \mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹, and J-OFURO 453×10³ km³ yr⁻¹. Note that HOAPS 1σ uncertainty is as large as the range among satellite-based E_ocean and more than 3 times the corresponding SD, again implying an overestimation of the HOAPS uncertainty range (see Sect. 4.2). The OAFlux 1σ uncertainty is of the same magnitude as the SD among satellite-based E_ocean, whereas the SEAFLUX uncertainty estimate is small in comparison. The small E_ocean found by HOAPS is partly due to data coverage, as data are only available over the ice-free ocean within 80^∘ of the Equator. A test with ERA5 data showed that E_ocean decreases by 5 % when the data are adapted to HOAPS coverage. Conversely assuming a 5 % increase for HOAPS yields 417×10³ km³ yr⁻¹. The same reasoning applies to the other satellite data sets with similar effects on E_ocean.

The spread in P_ocean is of the same magnitude as that found for E_ocean: HOAPS yields $\mathrm{335}\pm \mathrm{44}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹, and GPCP $\mathrm{384}\pm \mathrm{31}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹, assuming uncertainty ranges of 13 and 8 % for HOAPS and GPCP, respectively. ERA5 yields $\mathrm{426}\pm \mathrm{2}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹, which is significantly larger than either HOAPS or GPCP. From the GPCP-1DD data used in this study, we determine $P_{\text{land}}=\mathrm{116}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹, which is close to the estimates presented previously (which range from 111×10³ km³ yr⁻¹ (Oki and Kanae, 2006) to 117×10³ km³ yr⁻¹ (Rodell et al., 2015)). This is due to the fact that GPCP is used for all observation-based estimates. ERA5 P_land is somewhat higher than the observations (122×10³ km³ yr⁻¹), but the difference is not significant.

From the estimates of E_ocean and P_ocean it follows that for HOAPS $(E-P{)}_{\text{ocean}}=\mathrm{65}\pm \mathrm{106}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹, and for OAFlux $(E-P{)}_{\text{ocean}}=\mathrm{35}\pm \mathrm{48}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹. IFREMER yields $(E-P{)}_{\text{ocean}}=\mathrm{38}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹, J-OFURO $(E-P{)}_{\text{ocean}}=\mathrm{59}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹, and SEAFLUX $(E-P{)}_{\text{ocean}}=\mathrm{61}\pm \mathrm{31}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹.

The spread in R (40–51×10³ km³ yr⁻¹) is quite large. The total continental runoff from ERA5 is 41.2×10³ km³ yr⁻¹, which is slightly higher than the 39×10³ km³ yr⁻¹ found in ERA-Interim (Berrisford et al., 2011). Data from the Global Runoff Data Centre (GRDC; Wilkinson et al., 2014) yield an average of 41×10³ km³ yr⁻¹ with a standard deviation of 1.8×10³ km³ yr⁻¹ for 1987–2010. They are at the lower bound of the estimates by Clark et al. (2015), who find $\mathrm{44.2}\pm \mathrm{2.7}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹ for 1950–2008. The same study cites estimates by various authors that range from 25 to 50×10³ km³ yr⁻¹, with those based on freshwater fluxes representing the lower boundary (25–39×10³ km³ yr⁻¹). The long-term average runoff estimated from the GRUN (Global RUNoff; Ghiggi et al., 2019) data set is 38×10³ km³ yr⁻¹, consistent with the abovementioned range, albeit somewhat smaller than the best estimate by Clark et al. (2015). Note that GRUN runoff estimates are not independent of reanalysis data, as the machine-learning algorithm uses surface temperature and P data from the 20th Century Reanalysis (Compo et al., 2011; Ghiggi et al., 2019). Improvements in the quality of E−P estimates will aid the quantification of river runoff by providing an independent estimate of the total freshwater flux.

Where runoff is the net transport of (liquid) water from land to ocean, over-ocean VIMD (∇⋅(vq)_ocean) is, to a good approximation, the net amount of water vapor advected from ocean to land. Hence, $R=\mathrm{\nabla }\cdot (vq{)}_{\text{ocean}}=-\mathrm{\nabla }\cdot (vq{)}_{\text{land}}$. Whereas ERA5 estimates of ∇⋅(vq)_land are at the high end of the range of R mentioned above, at 31–33×10³ km³ yr⁻¹, ∇⋅(vq)_ocean is too small. As observed above, the fact that ERA5 VIMD is calculated in grid point space causes global total ∇⋅(vq) to be about 10×10³ km³ yr⁻¹, not 0. In addition, due to the tighter observational control over land, analysis increments may be larger over ocean than over land and may cause net ∇⋅(vq) to be close to net E−P over land, but less so over ocean (Paul Berrisford, personal communication, October 2020). There is another field in the ERA5 archive, the vertical integral of divergence of moisture flux (VIWVD, parameter ID p84.162), which is very similar to VIMD but is computed from hourly instantaneous reanalysis fields (and contains no contributions from liquid or solid water – but these can be neglected for our purposes). Globally VIWVD adds up to 0.9×10³ km³ yr⁻¹ (0.003 mm d⁻¹), a factor of 10 smaller than total VIMD. In addition, the agreement between over-ocean VIWVD and (E−P)_ocean is much better than that found for VIMD and (E−P)_ocean, and at 41.6×10³ and $-\mathrm{40.7}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹, respectively, over-ocean VIWVD and over-land VIWVD are also in agreement with other values in the five rightmost columns of Table 4. The estimates of net transport by Oki and Kanae (2006), Trenberth and Asrar (2014) and Rodell et al. (2015) are in agreement with ERA5 R and ∇⋅(vq)_land. The consistency between runoff and net transport seen in the last four rows of Table 4 is mainly by construction, as both are usually required (or defined) to be equal.

The five rightmost columns of Table 4 should, theoretically, all contain identical values (except for the sign). In practice, however, (E−P)_land ranges from −46 to $-\mathrm{40}\times {\mathrm{10}}^{\mathrm{3}}$ km³ yr⁻¹, (E−P)_ocean from 24 to 52×10³ km³ yr⁻¹, ∇⋅(vq) from 40 to 51×10³ km³ yr⁻¹ (disregarding the erroneous ERA5 ∇⋅(vq)_ocean value), and R from 40–51×10³ km³ yr⁻¹. Assuming the degree of consistency found among these values represents the reliability of the estimate, it is clear that E−P uncertainty is largest over ocean, and from the first two columns of Table 4 it follows that E and P contribute almost equally to that uncertainty.