Global hydro-climatological indicators and changes in the global hydrological cycle and rainfall patterns

There are few commonly used indicators that describe the state of Earth’s global hydrological cycle and here we propose three indicators to capture how an increased greenhouse effect influences the global hydrological cycle and the associated rainfall patterns. They are: i) the 24-hr global total rainfall, ii) the global surface area with daily precipitation, and iii) the global mean precipitation intensity. With a recent progress in both global satellite observations and reanalyses, we can now estimate the global rainfall surface area to provide new insights into how rainfall intensity changes over time. Based on the ERA5 reanalysis, we find that the global area of daily precipitation decreased from 43 to 41% of the global area between 1950 and 2020, whereas the total daily global rainfall increased from 1440 Gt to 1510 Gt per day. However, the estimated 24-hr global precipitation surface area varies when estimated from different reanalyses and the estimates are still uncertain. To further investigate historical variations in the precipitation surface area, we carried out a wavelet analysis of 24-hr precipitation from the ERA5 reanalysis that indicated how the rainfall patterns have changed over time. Our results suggest that individual precipitation systems over the globe have shrunk in terms of their spatial extent while becoming more intense throughout the period 1950–2020. Hence, the wavelet results are in line with an acceleration of the rate of the global hydrological cycle, combined with a diminishing global area of rainfall. the annual mean estimates. Brief ‘negative’ spikes can be seen for January 1st in 1950 and 1979 when the two sets of integrations started and did not calculate a full 24-hr cycle of precipitation. These results indicate that the global surface area with 24-hr precipitation has decreased between 1950 and 2020 from 0.43 to 0.41.


Introduction
The most common index for the state of Earth's climate and global warming has traditionally been the global mean temperature or the global sea level [1]. The global mean temperature has played a key role in determining the climate sensitivity and has been used as a yardstick for the Paris Agreement for limiting climate change to 1.5˚C or 2.0˚C warming from preindustrial times [2]. The global mean sea level represents an indicator that integrates the total effect of global warming, and functions like the mercury in a thermometer, as the ocean volume expands with higher temperature, in addition to a contribution from melting glaciers and ice sheets. One advantage of using the global mean sea level as a measure of the state of the climate system is that it is not biased by an irregular distribution of thermometers [3] nor the urban heat island effect [4]. Nevertheless, climate change is more than changing temperatures, and a strengthened greenhouse effect is also expected to change the hydrological cycle [1,[5][6][7][8] which can be described as a chain of processes connected in a closed-loop where H 2 O exists in various forms within the planetary system. The hydrological cycle is also closely connected to the flow of energy in the atmosphere through latent heat uptake (evaporation) and release associated with convection, atmospheric overturning, and condensation [9][10][11], and it produces typical rainfall patterns with wet and dry regions characterised by an unsteady appearance that may vary between the extreme states of droughts and floods depending on location. Furthermore, the surface temperature is affected by both cloudiness, soil moisture, and vegetation, all of which are influenced by prevailing rainfall patterns. Hence, the nature of the hydrological cycle is a key factor for shaping Earth's climate and has profound consequences for ecosystems and society. It is therefore important to summarise its state through a set of relevant indicators. Despite its profound consequence for Earth's climate, there is yet no widely used global indicator to describe the state of the global hydrological cycle. For instance, the global climate indicators of the World Meteorological Organisation (WMO) https://gcos.wmo.int/en/globalclimate-indicators include no measure of the hydrological cycle. In the IPCC's sixth assessment report (AR6) working group 1 (WG1), changes in the global hydrological cycle are represented by a set of indices, namely precipitation-evaporation differences over ocean and land (P-E), hydrological sensitivity (η = 2.1 − − 3.1%/˚C), total column water vapour, surface humidity (specific and relative), and global river runoff, all of which are discussed in greater detail in Chapter 8 [2]. However, Chapter 1 of AR6, which both sets the scene for the entire WGI assessment and presents the main concepts and methods, only discusses the CO 2 concentrations, the global mean temperature, and the global mean sea level, but no indicators for the state of the global hydrological cycle, such as the total global 24-hr rainfall, the global area with rainfall, or the mean precipitation intensity, hence confirming that global hydro-climatological indicators are not commonly used, just as the WMO list suggests. Furthermore, our ability to estimate such indicators has recently improved with new global data products becoming available, including both state-of-the-art high-resolution reanalyses like ERA5 [12,13] and remote sensing data from satellites. We need large volumes of daily global precipitation data to estimate the total area of rainfall on a daily basis, but once these indicators are estimated, they serve to aggregate the state of the hydrological cycle in terms of extremely small data volumes. There have already been some suggestions to include hydrological-based indicators such as the 'hydroclimatic intensity' [14] and the global precipitation area A p [6,15]. The latter is an attractive index since it represents a physical constraint connected to the mean precipitation intensity μ, which can be taken as the ratio of the total 24-hr rainfall amount on Earth P t for a random day to the global area of precipitation A p : μ = P t /A p . The value of P t is determined by the integrated global rate of evaporation E, as the global hydrological cycle involves a continuous and steady global flow of H 2 O through the climate system. In other words, what goes up must come down, and for a steady planetary system in equilibrium with a closed-loop of exchange of H 2 O between the surface and atmosphere, the evaporation is balanced by the precipitation according to P t = −E. A consequence of the closed circulation is that a smaller precipitation surface area implies more concentrated rain, with higher intensity, and less frequent rainfall over parts of the planet [15]. Furthermore, global evaporation is expected to increase with global warming [2], and hence the global surface area of daily rainfall A p and the total daily global rainfall amount P t may be considered as two key aspects expected to influence the statistics of heavy rainfall.
Here we present an analysis of three key indicators that aggregate information about the global hydrological cycle that includes the scale of rainfall patterns as well as their intensity: i) the 24-hr global precipitation surface area, ii) the total global 24-hr rainfall, and iii) the global mean precipitation intensity. We used a threshold of 1 mm/day to distinguish between dry and wet grid boxes. An analysis of the global 24-hr rainfall surface area based on the Tropical Rain Measurements Missions (TRMM) has indicated that the semi-global (50˚S-50˚N, representing 77% of Earth's surface) rainfall area diminished by 7% over the period 1998-2016 [15]. State-of-the-art reanalyses may also be used to estimate the complete global rainfall area, and here we present an updated analysis for the period 1950-2020 based on the ERA5 reanalysis [12]. The reanalyses provide the best information that we have on the atmosphere but come with several caveats. In a recent study by Bandhauer et al. [16], the ERA5 reanalysis was compared against observational gridded datasets for Europe that had been derived through statistical interpolation of rain-gauge observations. The comparison was carried out on a continental scale, using the daily gridded observational dataset E-OBS [17] as well as more regional high-resolution datasets over the three subregions: the Alps, the Carpathians, and Fennoscandia. The ERA5 reanalysis was found to reproduce the pattern of daily precipitation climate in all three sub-regions, and agree reasonably well with the E-OBS predictions. Bandhauer et al. [16] concluded that the ERA5 is consistent with independent observations at the mesoscale, but has some systematic biases in terms of overestimating the precipitation amounts in all regions, particularly because of the simulation of too many wet days. The bias was found to be more significant during summer.
To underscore the value of these global hydrological indicators and provide a detailed account of the conditions that they summarise, we conducted a more comprehensive study of the spatial characteristics of 24-hr precipitation using wavelet analysis. The global surface area with precipitation is the sum of the spatial extent of simultaneous precipitation generating systems at various spatial scales. Wavelets characterise the frequency and intensity of the precipitation features as a function of their spatial dimension, and hence can describe how their spatial extent has varied over time. Wavelet transforms are described in detail by Mallat [18], and wavelet decomposition is often used in spatial forecast verification to assess the skill of numerical model predictions on different spatial scales [19,20]. The multi-resolution analysis performed in this study was based on that proposed by Casati [20], but here we applied the wavelet transform directly to daily precipitation fields (without thresholding) as in Lussana et al. [21]. The daily gridded precipitation fields were decomposed onto spatial scale components by a 2-D Haar discrete wavelet transform, and their energy spectra were analysed. For each spatial scale, the wavelet energy is proportional to the number of precipitation features and their intensity, whereas the wavelet energy percentage provides information about how they are distributed across the scales, and hence provides information about changes in the scale structure of the precipitation fields. Further details about the analysis and data processing are explained in the Methods section in the S1 Appendix.

Results
Hydro-climatological indicators, such as the global rainfall area and global rainfall totals, are sensitive to the geographical distribution of persistent rainfall and the typical intensity of the rainfall. It is therefore instructive to start with maps showing how the mean and trends depend on location. Fig 1 presents maps of average wet-day mean precipitation m � , estimated from the ERA5 reanalysis, and its trend estimates over the 1950-2020 period (in terms of proportional change per decade) with only the statistically significant estimates being highlighted. We used the Students t-test at the 0.05 significance level to identify statistically significant trends. High values of m � (>10mm/day) are found along the equator which can be explained by the Inter-Tropical Convergence Zone (ITCZ). There are also moderately high values of m � in the region of the Gulf Stream extension and the Kuroshio in association with the storm track regions. The lowest values of m � (<2mm/day) can be found over the relatively cold maritime regions in the southeast Pacific and the southeast Atlantic, dominated by low stratocumulus clouds. According to the trend analysis for μ, there has been a general increase across most of the globe, particularly in the southern hemisphere and along the equator in the Pacific. A general increase in μ is also consistent with previous trend analysis based on rain-gauge data [22]. Our trend analysis indicated that 78% of the global area had an increase over the 1950-2020 period and 76% between 50˚S-50˚N. In this global summary, we counted all trend estimates and subjected them to a Walker's test that indicated that they also exhibited field significance on a global basis [23]. On a regional The wet-day mean precipitation μ is a statistic calculated for each year by selecting 'wet days' (days with more than 1 mm recorded precipitation) and estimating their average value, and is also referred to as 'mean precipitation intensity'. Both mean and trends were estimated from these annual statistics (The coastline data ETOPO1 from https://www.ngdc.noaa.gov/mgg/global/, doi:10.7289/V5C8276M).
https://doi.org/10.1371/journal.pclm.0000029.g001 basis, statistically significant decreases have taken place over parts of central Africa, the southeast Pacific, and west of India. An apparent strong decrease is seen over Sudan/Ethiopia, but it is a result of low mean values for m � in the denominator of the proportional estimates.  decrease in f w , and there has been a particularly pronounced reduction over the already dry regions in Africa. If we limit the analysis to 50˚S-50˚N, as in the TRMM data in [15], then 74% of the area exhibited a decrease in f w over the period 1950-2020. These results may also be connected to a reduction in the global surface area of rainfall.
Reanalyses make it possible to estimate the global sum of precipitation per day P t in terms of gigatonnes (Gt) and based on ERA5, P t has increased from 1440 Gt/day in 1950 to 1511 Gt/ day in 2020 (based on the linear trend fits shown as dashed lines in Fig 3). This trend is statistically significant at the 1% level, and consistent with an increase in the global rate of evaporation associated with higher global mean surface temperatures. The global mean surface air temperature estimated from ERA5 for 1950 was 286.8 K and 287.9 K for 2020 (13.7 and 14.7˚C respectively), implying an evaporation rate that increased by 5% per degree warming if we can assume an equilibrium state where E = −P t . This change is also consistent with increased latent heat flux near the surface and an accelerated turnaround rate in the global hydrological cycle.
The global daily precipitation surface area A p of ERA5 showed a statistically significant (at the 1% level) negative trend over the period 1950-2020 (Fig 4; from 43% of the global surface area in 1950 to 41% in 2020). Most of the reduction over time appears to be associated with a jump between 1985 and 1990, which is seen both in the fraction of the global surface area with rainfall and the rainfall area between 50˚S-50˚N (Fig 4). One hypothesis is that the jump in A p estimated from the ERA5 was an artifact associated with the assimilation of different satellite instruments. Similar estimates of the precipitation surface area from a range of reanalyses (NOAA 20CRv3 [24], ERA20C [25] and NCEP1 [26]) showed a different mean state and Daily and annual total global rainfall in Gt/day estimated from the ERA5 reanalysis. The blue curve represents the rainfall area between 50˚S and 50˚N. The thick curves show the annual mean P t estimated from daily estimates. Dashed lines are best-fit linear trend fits from an ordinary linear regression on the annual mean estimates. Brief 'negative' spikes can be seen for January 1st in 1950 and 1979 when the two sets of integrations started and did not calculate a full 24-hr cycle of precipitation. These results indicate a long-term increase in the total global precipitation amount from 1440 Gt/day in 1950 to 1511 Gt/day in 2020.
https://doi.org/10.1371/journal.pclm.0000029.g003 evolution: both NOAA and NCEP1 indicated statistically significant decreases while the rainfall area of ERA20C suggested little change over time (Fig A in S1 Appendix-ERA20C showed a tiny but significant decrease in rainfall between 50˚S and 50˚N, but an increase in global 24-hr rainfall surface area). The reanalyses differed in terms of spatial resolutions and covered different time periods, which may have influenced their representation of the rainfall patterns, their global structure, and trends. The NCEP1 reanalysis indicated a rapidly diminishing precipitation surface area from a fraction of 0.45 to 0.42 over the period 1961-2020, while the ERA5 showed a similar evolution in the 1990s. For ERA20C, the trend analysis revealed small but statistically significant changes (p < 0.001) over the period 1961-2015: an increase in the global rainfall area (from 0.441 to 0.444), but a decrease between 50˚S-50˚N (from 0.453 to 0.450). The NOAA 20CRv3 data showed a strong decrease throughout the whole reanalysis period, however, the high area fraction of daily rainfall derived from it was unrealistic (0.60-0.54 in the period 1961-2015; see Fig A in S1 Appendix). The differences between the reanalyses may suggest that the precipitation surface area is not closely constrained by the assimilation of observations, as the NOAA 20CRv3 and ERA20C reanalyses didn't involve assimilation of satellite retrievals. Furthermore, a comparison between the aggregated results derived from the TRMM and ERAINT used by Benestad [15] with the more recent ERA5 reanalysis revealed a difference in the rainfall surface area estimated from the two data sources. The analysis of the ERA5 reanalysis yielded an estimate for A p between 50˚S-50˚N of the order of 40%, whereas the TRMM data only gave an area of about 24% (see Figs A and B in S1 Appendix). The discrepancy in A p between TRMM and the reanalyses needs to be explained, and one potential reason may be the difference in spatial resolutions.
The observed increase in the global total mass of precipitation over a reduced global surface area has resulted in greater global mean precipitation intensities (Fig 5), which from the ERA5 data was estimated to be 6.1 mm/day in 1950 and 6.8 mm/day in 2020, with a trend that was statistically significant at the 1% level. Hence, increases in extreme 24-hr precipitation amounts can be explained in terms of both increased evaporation (more H 2 O into the atmosphere and increased circulation within the global hydrological cycle) and a reduced global rainfall surface area in ERA5, which may be due to dynamical changes that reduce the spatial extent of the precipitation generating systems and/or changes in the observational input.
The global hydro-climatological indicators reflect the sum of changes taking place on regional and local scales which were investigated with a wavelet multi-resolution analysis. There has been a shift in the wavelet energy of daily precipitation towards higher values for all spatial scales, indicating a general increase in rainfall. The change in the spatial characteristics of the hydro-climatological systems can be seen in Fig 6 which compares the statistics of the multi-resolution analysis of global daily precipitation fields over the two normal periods 1961-1990 (blue) and 1991-2020 (red and pink), where the (approximate) spatial scale for the decomposition level is shown on the abscissa while the wavelet energy is shown on the ordinate. Each spatial scale can be interpreted as a measure of the size of precipitation features. The results of the multi-resolution analysis are also summarised in Tables 1 and 2. The spatial scales smaller than 0.2 degrees may have been influenced by our choice of using bilinear interpolation to regrid ERA5 fields onto the dyadic domain, and for this reason, they are not shown in Tables 1 and 2. The first column of Table 1 shows a conversion table between degrees and km, where the kilometers are approximated to the nearest ten. For each spatial scale, the distribution of wavelet energy in the two normal periods is shown through boxplots (Fig 6). The values of the medians of the boxplots are presented in Table 1 as q 50 ðEn 2 l Þ over 1961-1990 and 1991-2020. The proportional contributions of the spatial scales are presented in Table 2 in terms of percentage, and the meso-β scale between 20 km and 200 km represents roughly onethird of the total energy of the wavelet coefficients, while atmospheric phenomena on the macroscale represent the remaining two thirds. The peak in the wavelet energy intensities was found for the scale of 4 degrees (440 km), and the three scales of 2, 4, and 8 degrees (220 km to 890 km) alone contain about half of the total wavelet energy.
One remarkable result was that the median values for spatial scales between 0.5 to 8 degrees during 1991-2020 were higher than the 75th percentile of the same spatial scales found for 1961-1990, thus indicating a significant shift in intense precipitation features at these scales.  distribution. On the other hand, outliers indicated a small number of individual days with exceptionally high or low energies during both normal periods. For instance, the single day with the highest wavelet energy value belongs to the normal period  and was reached at the spatial scale of 4 degrees. In Table 1, the ratio between q 50 ðEn 2 l Þ at 1991-2020 and 1961-1990 shows the relative changes in wavelet energy at each spatial scale. The relative increment was larger for the smaller spatial scales, and the signs of the relative changes in Table 2 show that there has been an energy transfer between spatial scales, with the energy of the wavelet coefficients transferring from synoptic to mesoscale. In other words, the wavelet analysis applied to the ERA5 reanalysis indicated that typical mesoscale atmospheric phenomena, such as squall lines or thunderstorm groups, have become more important in determining the total wavelet energy balance of daily precipitation. The same multi-resolution analysis was performed separately for three regions over the globe: the northern and southern temperate zones (TN and TS, respectively) and the tropics (TR). This analysis suggested that the variations Fig 6 were not occurring uniformly over Earth's surface, rather the change was much more pronounced in the tropics than in the temperate zones (see Fig C and Table A in S1 Appendix), as the trends presented in Figs 1-2 also seem to suggest.
Time series of wavelet energies at different spatial scales between 1950 to 2020 indicated a gradual shift in the wavelet energy distribution across the spatial scales over time (Figs A and Table 1. Summary of the comparison of the multi-resolution decomposition wavelet energies over the entire globe between 1961-1990 and 1991-2020. q 50 (En 2 l ) is the median of the energies at spatial scale l. The relative change is the ratio between the energy medians in 1991-2020 and in 1961-1990, multiplied by 100. The linear trend is the angular coefficient of the best-fitting line of the daily energies in the period 1991-2020. Note that only spatial scales greater than or equal to 0. 25  D-I in S1 Appendix). The shift took place at an approximately constant rate after 1985, however, the time series were flat before then. From around 1985 and onward, the wavelet energies at all scales began to increase, each spatial scale with a different growth rate which is presented as linear trends over the period 1991-2020 in Table 1, and in terms of daily wavelet energy fractions in Table 2. At the Meso-β scale and the lower part of the synoptic scale (up to 440 km), the linear trends of the daily wavelet energy fractions were positive (ranging between 1 and 14%; Table 2), thus remarking again that this part of the energy spectrum has become more important in the total wavelet energy balance of daily precipitation. Similar time series for the wavelet energies in the tropics, northern and southern temperate zones indicated similar shifts in the wavelet energies in the tropics as for the global domain (Fig E-G in S1 Appendix). However, the wavelet energies in the tropics decreased from 1950 to 1985, and hence differed from those for the entire globe, but after 1985 they also began to increase over time and with a faster growth rate for smaller spatial scales. In the northern temperate zone, there has been a gradual increase in the precipitation wavelet energy for all spatial scales, albeit less pronounced than in the tropics. The evolution of the wavelet energy time series in the southern temperate zone was partly similar and partly different from other areas: As for the other two regions, the wavelet energies in the southern temperate zone increased over all spatial scales, but with the clear distinct feature that the increase was gradual over the whole 1950-2020 period. Another difference was that the wavelet energies of the largest spatial scales grew faster than those of the smallest spatial scales, and for this reason, the daily percentage of wavelet energy also increased for the larger spatial scale (Fig G in S1 Appendix). The growth rates in the tropics were at least one order of magnitude greater than in the temperate zones (Fig E in S1 Appendix), and hence underline the importance of the changes in the rainfall patterns taking place at low latitudes. We applied a linear regression to the annual mean energy of each global wavelet component and the global mean temperature, also estimated from ERA5. The results from this regression analysis indicated that the trend over time in wavelet energy closely followed the global mean temperature (Fig 7), albeit with some deviations around 1980-1990. The strikingly similar evolution in both the global mean temperature and wavelet components 4-8 and 12 may suggest that the global mean temperature influenced the spatial extent of precipitation features, hence providing an empirical connection between the increased greenhouse effect and changes in the global hydrological cycle. It is also in line with previously reported linear dependency between the global mean temperature and the 50˚S-50˚N daily precipitation area [15]. Such a dependency on the global mean temperature may enable downscaling of information about storm systems in a new way.
Fig 8 provides a case study illustrating how the shift in the wavelet energy distributions takes place in reality for the ongoing changing statistics of daily precipitation. Panel b shows the daily precipitation simulated by ERA5 in the Gulf of Mexico on the 28th of August 2005 in connection with the transit of Hurricane Katrina at its peak intensity over that region. Hurricane Katrina was a Category 5 Atlantic hurricane that caused enormous damages to the regions it crossed. The upper part of each panel shows the map of daily precipitation, while the bottom part shows the precipitation intensity as a function of the longitude to get a horizontal "profile" of precipitation intensity. The multi-resolution analysis enabled us to obtain relative changes between the wavelet energies in the normal periods 1991-2020 and 1961-1990, as shown in Tables 1 and 2. The ratio between the medians of wavelet energies corresponding to the same spatial scales as those reported in Table 1 over the two generic time periods A and B (i.e. median of the energy on A divided by median on B) expresses how the variability of the wavelet coefficients changes, scale by scale. If this ratio is higher than one for a specific spatial scale, then it is more likely to have larger values of the wavelet coefficients on that scale during period A than during period B. If this ratio is less than 1, however, then the it is more likely to have smaller values. Therefore, the relationships between the medians of the wavelet energies can help us to estimate what an event would have been like if it happened in a different period than the one in which it actually happened. For the analysis of Hurricane Katrina, the discrete wavelet transformation was applied only to the precipitation field over the entire globe on the 28th of August 2005. It was then possible to re-scale the wavelet coefficients using the square root of the ratios between q 50 (En l ) 1961-1990 and 1991-2020 (i.e. the inverse of the relative deviations reported in Table 1) as weights. For instance, with reference to Table 1, the wavelet coefficients for scale l = 2 degrees have been multiplied by 0.89, that is the square root of 4.55/5.72, to transform them from the actual ones observed on the 28th of August 2005 to those of a hypothetical twin of Hurricane Katrina in the 1961-1990 climate. The overall effect was to "deflate" the wavelet energies in a way that was consistent with the observed average variations in the medians of the distributions over the two normal periods. The precipitation field for Hurricane Katrina on the 28th of August 2005 adjusted to the 1961-1990 climate is shown in Panel a. The growing trends in the wavelet energies shown in Fig 7 from 1991-2020 and in S1 Appendix was in this case used to extrapolate the future wavelet energies for all days in the next normal period 2021-2050, assuming a constant growth rate in the wavelet energy that continues unaltered over the whole period (i.e. an increasing monotonic linear trend of the wavelet energies). For each wavelet component, the slope and intercept coefficients of the best-fit line for the daily wavelet energies over the period 1991-2020 were computed. Then, the coefficients were used to linearly extrapolate daily wavelet energies over the period 2021-2050 and compute corresponding statistics for 2021-2050 as those reported in Table 1. Once those scale-dependent future energies were obtained, their relative changes in terms of their medians with respect to 1991-2020 were computed as we did in Table 1 for other periods. At this point, we re-scaled Hurricane Katrina's wavelet coefficients as we had done above for the 1961-1990 climate, but this time we got a reconstruction valid for the period 2021-2050. In this example, the wavelet coefficients for scale l = 2 degrees was inflated since they were multiplied by a factor of 1.14. Such a case is shown in Panel c.

Discussion
Here, we have presented a 2-D wavelet analysis of the 24-hr precipitation from the ERA5 reanalysis to distill new details about the nature of weather systems producing precipitation that may affect the global hydrological indicators. It is interesting to relate these wavelets to the daily global surface area of precipitation greater than 1 mm, as estimates based on ERA5 illustrate how the rainfall patterns over the period 1950-2020 may have become more intense and shifted towards smaller spatial structures. These changes are consistent with the proposed aggregated summaries of the hydrological cycle: the total global 24-hr precipitation mass, the global precipitation surface area, and the global mean precipitation intensity. It is possible that a reduction in the global area of precipitation is linked to a deceleration in the motion of the rain-producing systems [27], or that the drop in A p between 1980 and 1990 may be related to aerosols and a global dimming [28] or an artifact of changes in the assimilated data (being timed with discrepancies seen in Fig 7), but it could also be due to a change in the cloud structure in line with the multi-resolution analysis. While none of these explanations exclude each other, we note that the latter is consistent with a simple conceptual representation of an increased greenhouse effect entangled with the global hydrological cycle, with increased convective and latent heat transport in an atmosphere that is more opaque to infrared radiation [11]. The results of the 2-D wavelet analysis also suggest that the rain-producing systems have reduced in spatial extent in the tropics, but probably not at higher latitudes. Hence, these wavelets are not necessarily inconsistent with results from a recent analysis of regional climate models from the Euro-CORDEX ensemble, which suggested that large precipitation systems over Europe have become more frequent and even larger with global warming, while precipitation systems of lesser extent will be reduced in numbers [29]. One difference between those findings and ours is that they only involved events corresponding to the 90th and 99th percentiles, whereas we analysed the area of common rain events with a threshold of 1 mm/day.
Our results also leave us with some new questions such as whether differences between hemispheres were due to different coverage of observation or perhaps different land areas. Even if the trends found were artifacts of the reanalysis setup, it is important to document them, and if they give a false impression, then this would be the first step in making further progress. One issue is that significant variations in the observational network feeding a fixed model with data for assimilation to generate all ERA5 precipitation fields, such as the massive increase in the amount of satellite data assimilated, are likely to have a measurable impact on the quality of the simulations. Substantial variations in A p from different reanalyses may be a demonstration of the effect of data assimilation. In this case, we assume that ERA5 provides a better representation of Earth's climate system as it's more recent and advanced than the others, although their differences provide some indication of the level of uncertainty. The amount of observations used in the ERA5 data assimilation cycle grows over time, especially after the beginning of the satellite era which started approximately in 1980. On the other hand, ERA5 didn't assimilate precipitation observations directly, but only observations on the vertical profile of humidity, and the volume of such observations is changing over time. There was one local exception over North-America after around 2008 according to [12]: "ERA5 assimilates the NCEP stage IV quantitative precipitation estimates produced over the USA by combining precipitation estimates from the NEXRAD with gauge measurements". Another potential caveat is that ERA5 was provided in two parts: ERA5 back extension (preliminary version) from 1950-1978 and ERA5 1979-2020; and the beginning of the increase in the running mean of the wavelet energies overlaps in part with the type change of the ERA5 dataset. Furthermore, the dates 1950-01-01 and 1979-01-01 ERA5 do not include precipitation for the complete day since six hours spin-up time was excluded in the reanalyses. This minor inconsistency does not affect the long-term trends but is visible as brief 'negative spikes' in Figs 3 and 4. Nevertheless, one limitation with hydro-climatological indicators based on the ERA5 reanalysis is that the estimates can potentially be connected to the increase in the number of observations assimilated (more and more satellite data). The finer the spatial resolution of the overall observational network, the "sharper" the precipitation events can be simulated. When more data are assimilated, the effective resolution of the analysed field becomes finer, and the final result is that the precipitation fields are less blurry. The availability of observations for assimilation can also explain why ERA5 reanalysis has large uncertainties over Africa because there are few observations available for the assimilation of the atmospheric models. The coincidence between the beginning of the gradual increase in wavelet energies and the beginning of the satellite era, however, gives rise to the suspicion that the climatic signal was distorted by variations in the observational network.
Another caveat connected to global aggregates of rainfall is that neither the reanalyses nor the satellite-borne remote sensing necessarily provide an accurate representation of precipitation at spatial scales smaller than the mesoscale [16]. These data sources are nevertheless the best we have, and they do offer useful information about the rainfall patterns. It is also important to bear in mind that other measurements also have limitations, e.g. the global mean temperature [3,30]. Albeit imperfect, a more widespread use of such indicators may prompt further efforts to improve their estimation over time. One nagging question is why the estimates of A p differ between TRMM [31,32] and ERA5 (Figs A and B in S1 Appendix). One plausible explanation is that TRMM measures instantaneous conditions during each satellite overpass, sampling the same spot only a few times per day (approximately 3hr intervals), while ERA5 accumulates precipitation over 24 hrs. Since the rain generating systems often last minutes to hours, we should expect a larger area when we accumulate over 24 hrs as opposed to the daily sum of snapshots approximately every 3 hr depending on the satellites' times of overpass. The two datasets also have different spatial resolutions, where the TRMM data consists of processed retrievals with a 0.25-degree resolution whereas ERA5 has a T639 grid which implies a spatial resolution of 0.28 degrees. Another difference is that ERA5 is a product of model analysis calculated for each grid box, while TRMM is a gridded product from satellite retrievals. However, neither the reanalyses nor the satellite-borne remote sensing may necessarily provide an accurate representation of precipitation at spatial scales smaller than the mesoscale [16]. Orlanski [33] proposed a subdivision of scales for atmospheric processes based on their horizontal scale length, which has been reviewed by Thunis and Bornstein [34]. The term synoptic-scale or macroscale is used for atmospheric scales greater than 200 km and includes atmospheric phenomena such as synoptic cyclones, fronts, and hurricanes. The scales from 200 m to 200 km are included in the mesoscale. The choice of a specific aggregation period does influence the spatial scales that may be studied adequately without distorting the characteristics of the phenomena under examination. For daily aggregation, the lower bound of the resolved spatial scales is approximately 20 km and coincides with that part of the mesoscale that is named meso-β (i.e. from 20 km to 200 km). Phenomena within the meso-β comprise thunderstorm groups. The minimum lifetime of daily precipitation phenomena that we could study lies somewhere between a few hours and 1 day, and this limits the minimum spatial scale to the beginning of the Meso-β scale and approximately 20 km.
The multi-resolution analysis showed a shift in the distribution of daily wavelet energies between the two normal periods 1961-1990 and 1991-2020 towards higher energies for the most recent normal period, and we uncovered clues suggesting that the detected climate change in our results were indeed real. The most important among these clues was the fact that changes were seen at large spatial scales, such as precipitation patterns with length scales of 8 or 16 degrees, which should be well resolved already by a low-density observational network. The impact of an increase in the density of the observational network should be more evident for the smaller spatial scales, and almost non-existent for spatial scales larger than some degrees because the background atmospheric model used in the re-analyses for those scales harmonises the assimilated observations. Besides, our analysis on the southern temperate zone showed that the increase in precipitation wavelet energies took place over the whole time period considered, and the larger spatial scales have become more important over time than the smaller spatial scales. We also applied the same multi-resolution analysis to the ERA20C reanalysis [25] (see Figs H-I in S1 Appendix), which had a horizontal resolution of approximately 125 km and where the only observations assimilated included surface and mean sea level pressures and surface marine winds. The results from this reanalysis could not be affected by satellite retrievals because only in-situ surface observations had been assimilated, which also made this reanalysis an ideal dataset for assessing the influence of assimilated observations on the multi-resolution analysis results. The wavelet energies increased constantly over time for all spatial scales, however, the growth rates were different and depended on the spatial scale, and the energies of the smaller spatial scales grew faster than those of the larger spatial scales. The example shown in Fig 8 is just one way that allows us to show the impact of the variations in the wavelet coefficients on the precipitation fields in both graphic and quantitative terms. It is also possible to use the scaling and offset obtained from a regression analysis against the global mean temperature as shown in Fig 7 to 'downscale' future rain patterns. In addition, the close match between the temporal evolution in the global mean temperature and several wavelet components is also providing further support of veracity. Nevertheless, it is still uncertain whether these changes actually took place on Earth during 1950-2020, and there is a need for future efforts to elucidate the reasons for different results obtained between reanalyses as well as satellites.

Conclusions
We have proposed three global hydro-climatological indicators to provide a summary of the state of Earth's global hydrological cycle: (i) the global 24-hr total precipitation, (ii) the global surface area of 24-hr precipitation, and (iii) the global mean precipitation intensity. They provide an aggregated indicator of shifts in the rainfall patterns, while traditional trend analysis of wet-day frequency and wet-day mean precipitation estimated from the ERA5 reanalysis indicates changes in the precipitation statistics over the period 1950-2020 that can explain a decrease in the surface area of 24-hr precipitation and an increase in the intensity of the rainfall. Furthermore, a multi-resolution 2-D wavelet analysis applied to the same data has revealed changes in the intensity and frequency of precipitation systems, with different rates for different spatial scales (e.g. more dramatic for meso-β scales in the tropics) that appear to be correlated with the global mean temperature. Moreover, we found indications of a combination of an increasing trend in rainfall intensities and a shift in the presence of rainfall patterns from large to smaller spatial scales in the tropics, and (to a lesser extent) from small scales to large scales in the Southern Hemisphere extra-tropics. However, there was one important caveat, as these findings were subject to considerate uncertainties: there were substantial variations in estimated global precipitation surface area between different reanalyses.

Data & methods
The ERA5 reanalysis was produced in two chunks: 1950-1978 and 1979-2020, where the latter included assimilation of satellite retrievals. The ERA5 reanalysis hence may not be homogeneous due to the introduction of new and more advanced instruments over time [35]. This caveat should be borne in mind, while ERA5 is considered the best data we have for trying to quantify the global status on the hydrological cycle. More details about the methods are provided in S1 Appendix.
The ERA5 reanalysis was downloaded on a grid with a spacing of 0.25 degrees in both eastern and northern directions which approximately corresponds to 28 km if we assume the approximate metric equivalent for 1 degree to be 111 km. Spatial scales smaller than 28 km in the multi-resolution analysis were the result of regridding to higher resolution, and as expected did not exhibit statistically (not physically) significant signal.

Mathematical expressions
The global total rainfall P t (t) for any day t was estimated from the sum of the 24-hr precipitation over each pixel multiplied by the pixel area: The sum was taken over all n pixels/grid-boxes. We used a CDO command https://code. mpimet.mpg.de/projects/cdo to calculate the global average and the total global precipitation was taken as the product between Earth's surface area and the global averaged precipitation (the details are provided in the S1 Appendix). The global surface rainfall area A(t) was estimated in terms of the sum over pixel areas where 24-hr rainfall x i (t) exceeds the threshold of x 0 = 1mm/day for any day t: where Hð�Þ is the Heaviside function. Again, we used a CDO command that estimated the global mean and the estimate returned the fraction of Earth's area with 24-hr rainfall on day t (see S1 Appendix for the details), i.e. the fractional global surface area A(t)/A e (which also can be expressed in terms of percentages). The mean intensity was estimated according to

Trend estimates
The trend analysis presented here (e.g. Figs 1-5) involved an ordinary linear regression between annual aggregates (μ or f w ) and the year (τ) according tô whereŷ was a best-fit linear trend estimate. We also defined statistically significant trends as those where the estimate of β was different to zero with a statistical significance at the 5%-level (which implies a low estimate for the standard error δβ).

Wavelet analysis
The wavelet analysis performed in this study was based on the scale-separation verification approach proposed by Casati [20], with the sole difference that here we apply the wavelet decomposition directly on the precipitation fields without thresholding. In what follows, we outline the method and in S1 Appendix we provide the markdown of the R function used to perform the calculation. Wavelet transforms [18] are similar to Fourier transforms as they enable the decomposition of a signal (e.g. the precipitation field) into components with different frequencies (i.e. different scales). Wavelets are however locally defined, and hence result in more suitable transforms than Fourier for representing episodic and spatially discontinuous fields, characterised by onand-off features, such as precipitation. In our approach, we use 2-Dimensional (2D) Haar wavelets, which can be visualised as "square" sine functions.
The precipitation fields P are decomposed into the sum of components on different spatial scales by using a 2D Haar discrete wavelet transform where the mother wavelet components MW j on the scales j = 1, . . ., J are fields of effective resolution of 2 j−1 × 2 j−1 grid-points, and the largest scale FW J (named the "father wavelet" component) is the spatial average of the precipitation field over the wavelet dyadic domain of 2 J × 2 J grid-points. In order to simplify our notation, we index the scale components with r = 0, . . ., R = J, which is associated with their effective resolution of 2 r × 2 r grid-points. To perform the wavelet transform, the reanalyses precipitation fields were interpolated to a dyadic grid of 2048 × 2048 grid-points. The energy spectra of the precipitation fields is then evaluated by computing the spatial average of the squared values of all grid-points, for each separate scale component P r , r = 0, . . ., R, as: The energy En 2 is proportional to the number of precipitation features and their intensity, for each scale. The comparison of present versus future energy spectra, therefore, informs on the changes in precipitation occurrence and intensity, as a function of the scale of the precipitation features.
Due to the linearity of the wavelet transforms and the orthogonality of discrete wavelets, Casati [20] has shown that the energy of the un-decomposed (original) precipitation field P (which is the spatial average of the precipitation grid-point squared values) is equal to the sum of the energies of its scale components The energy percentage that each scale component P r , r = 0, . . ., R contribute to the total energy is then evaluated as En 2 ðP r Þ ¼ En 2 ðP r Þ=En 2 ðPÞ: ð8Þ The energy percentages inform how the precipitation features (occurrence and intensity) are distributed across the scales. Therefore, a comparison of present versus future energy percentages explains the evolution of the precipitation fields scale structure.

S1 Appendix. PDF-document with Figs A-I, Table A and output from R-markdown script.
(PDF) S1 File. Data file containing aggregated data used for the analysis. The raw data are also open access and access to those are explained in S1 Appendix, and the file also includes the Rcode (R-markdown) behind the analysis. It is also available through FigShare.com https:// figshare.com/ndownloader/files/31723382 [36]. (GZ)