Analysis of Water Yield Changes from 1981 to 2018 Using an Improved Mann-Kendall Test

Abstract

Water yield (WY) refers to the difference between precipitation and evapotranspiration (ET), which is vital for available terrestrial water. Climate change has led to significant changes in precipitation and evapotranspiration on a global scale, which will affect the global WY. Nevertheless, how terrestrial WY has changed during the past few decades and which factors dominated the WY changes are not fully understood. In this study, based on climate reanalysis and remote sensing data, the spatial and temporal patterns of terrestrial WY were revisited from 1981 to 2018 globally using an improved Mann-Kendall trend test method with a permutation test. The response patterns of WY to precipitation and ET are also investigated. The results show that the global multi-year mean WY is 297.4 mm/a. Based on the traditional Mann-Kendall trend test, terrestrial WY showed a significant (p < 0.05) increase of 5.72% of the total valid grid cells, while it showed a significant decrease of 7.68% of those. After correction using the calibration method, the significantly increasing and decreasing areas are reduced by 10.52% and 10.58% of them, respectively. After the correction, the confirmed increase and decrease in WY are mainly located in Africa, eastern North America and Siberia, and parts of Asia and Oceania, respectively. The dominant factor for increasing WY is precipitation, while that for decreasing WY was the combined effect of precipitation and evapotranspiration. The achievements of this study are beneficial for improving the understanding of WY in response to hydrological variables in the context of climate change.

1. Introduction

Water yield (WY) is the difference between precipitation and evapotranspiration (ET) and is an important indicator of regional water resource availability. Water forms the foundation of terrestrial ecosystems and is both the basic element of life and the site for a range of life activities [1]. In addition, water is a clean energy source [2,3,4,5]. Both greening and hydropower generation are restricted by water resources. Excessive changes in the availability of water resources may affect sustainable development in the local area. Thus, the WY changes in response to climate change caused by human activities should be seriously taken into account in this century for both humans and the biosphere [6].

According to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change [7], current global climate change is unprecedented. Climate change causes significant changes in the spatiotemporal distribution and circulation of precipitation and ET, which affect WY [8]. On the one hand, climate change exacerbates the global imbalance in precipitation patterns. Temporally, precipitation variability will increase with global warming on all time scales [9]. Spatially, precipitation variability is predicted to increase in the global humid zone (mainly in the tropical oceans, most monsoonal regions, and mid-to-high latitude regions). That is, the overall spatial pattern of precipitation changes shows a trend of “Dry gets drier, wet gets wetter” [10,11]. On the other hand, as an important part of the land surface hydrological cycle, ET, which comprises transpiration from plants and evaporation from the ground [12], is regulated by both climate change and land cover change. In general, meteorological factors were the dominant factors affecting ET; temperature, radiation, wind speed, and precipitation all have important effects on the interannual and seasonal variations of ET [13,14,15]. Further, vegetation cover also plays an important role in regulating ET [16]. Based on paired watershed experiments in small-scale watersheds, researchers have noted that vegetation restoration increases ET and canopy interception and, therefore, reduces runoff [17]. Conversely, other studies have shown that increased vegetation cover promotes the large-scale transport of water vapor and regional-scale precipitation; thus, vegetation restoration positively influences the hydrological cycle [18,19]. Considering the crucial role of WY and the complex effects of precipitation and ET, how terrestrial WY has changed during the past few decades urgently needs to be recognized for its prediction.

Previous studies have investigated the variations in WY in local areas [20,21]. However, few systematic analyses of global precipitation and ET data have been carried out [22,23]. This is not conducive to identifying and monitoring global WY changes under a unified standard framework. Linear regression is a commonly used method for estimating the parameters of long-series hydraulic data [24]. Compared to linear regression, the Mann-Kendall (MK) method does not need to meet a certain distribution to complete the statistical test. In addition, the structural system is less disturbed by outliers. Thus, it is more applicable and feasible than other methods that use a normal distribution test [25,26,27,28]. However, recent studies have found that the MK test has the problem of multiple hypotheses testing. This might lead to false positives and false interpretations of spatial patterns of change under the effect of spatial autocorrelation mechanisms [29,30]. In other words, ignoring the multiple hypothesis testing problems may overestimate the significance of the changes. According to a recent work by Cortés et al. [31], the multiple hypothesis test problem can be corrected using a cluster-based replacement test to clarify the spatiotemporal distribution of land surface greening. This proposed method provides an effective way to eliminate misjudged regions in the detection of long-term trends at regional and global scales.

Therefore, this study aimed to revisit the spatiotemporal patterns of terrestrial WY worldwide from 1981 to 2018. The WY was calculated using precipitation and ET based on long-term climate reanalysis data and remote sensing data. Significant changes in WY were identified using both the MK test and the method by Cortés et al. [31]. Furthermore, the study explored the response pattern of WY to precipitation and ET and the dominant factors of WY change under different climate types and vegetation distribution classification systems. This study will help to improve the understanding of WY in response to the hydrological variables in the context of climate change.

2. Materials and Methods

2.1. Precipitation Dataset

The precipitation data for this study were derived from the Climate Research Unit Time Series (CRU TS) v4.04 at the University of East Anglia (http://www.crudata.uea.ac.uk) (accessed on 20 October 2020). The dataset was derived by interpolating the monthly climate anomalies from extensive networks of weather station observations [32,33]. It includes several climatic variables at a 0.5° spatial resolution over all land domains of the world from 1901 to 2019, except Antarctica. In this study, monthly precipitation data from 1981 to 2018 were selected to calculate the annual precipitation in mm/a. This dataset has been widely used in global and regional climate change studies because of its complete coverage, high data resolution, and long-timescale [34,35,36].

2.2. Evapotranspiration Dataset

The land surface ET data used in this study were derived from the Global Land Evaporation Amsterdam Model (GLEAM), a global land evaporation product (http://www.gleam.eu) (accessed on 4 November 2019). GLEAM is a set of algorithms that separately estimate the different components of terrestrial ET based on satellite observations as follows: transpiration, interception loss, bare soil evaporation, snow sublimation, and open-water evaporation. Central to this methodology is the use of the Priestley and Taylor evaporation models [37,38]. In this study, the annual ET data (unit: mm/a) of the GLEAM v3.3a version with a spatial resolution of 0.25° × 0.25° from 1981 to 2018 were selected. To match the spatial resolution of the precipitation data, the ET data were aggregated to a 0.5° × 0.5° resolution.

2.3. Climate-Plant Types

In this study, each valid pixel was assigned to a climate-plant type, which was defined by both the Köppen climate classification and the Moderate-resolution Imaging Spectroradiometer International Geosphere Biosphere Program (IGBP) land cover type. The Köppen-Geiger (KG) climate classification was first developed in the late 19th century [39] and updated by Rudolf Geiger. The KG system classifies the climate into five main classes (Table 1). The classification was based on the threshold values and seasonality of monthly air temperature and precipitation. It is widely used and has been established as a multidisciplinary standard for characterizing the climate of a region [40].

The MODIS land cover type dataset (MCD12Q1) from 2001 at a spatial resolution of 0.5° was used to identify the land cover in this study. The IGBP classification system used in MCD12Q1 is an internationally common land use and land cover classification system. It defines 17 categories of land cover types, including 11 categories of natural vegetation, three categories of land use and land mosaic, and three categories of unvegetated land, ranging from 180°W to 180°E and 64°S to 84°N [41]. The IGBP land cover classification system reflects the characteristics of the physiological parameters of the land surface and the dominant vegetation status of the land cover. This is typical of a land cover classification system [42]. In this study, water regions were excluded.

Each of the climate-plant types was composed of climate and land cover types. For example, tropical evergreen coniferous forest was referred to as A-ENF, where A indicates a tropical climate and ENF indicates an evergreen coniferous forest (Table 1).

Table 1. Descriptions of the climate and plant types in this study are based on the Köppen-Geiger climate classifications and the MODIS IGBP land cover types, respectively. The details regarding the codes of Köppen-Geiger climate classifications are provided by Beck et al. [43].

Climate Types	A (Tropical)	B (Arid)	C (Temperate)
	D (Cold)	E (Polar)
Land Cover Types	Evergreen Needleleaf Forest (ENF)	Evergreen Broadleaf Forest (EBF)	Deciduous Needleleaf Forest (DNF)
	Deciduous Broadleaf Forest (DBF)	Mixed Forest (MIF)	Open Shrublands (OSH)
	Woody Savannas, Savannas (WSA)	Grasslands (GRA)	Croplands (CRO)
	Cropland and Natural Vegetation Mosaic (CNV)	Snow and Ice (SNI)	Barren or Sparsely Vegetated (BSV)

Table 1. Descriptions of the climate and plant types in this study are based on the Köppen-Geiger climate classifications and the MODIS IGBP land cover types, respectively. The details regarding the codes of Köppen-Geiger climate classifications are provided by Beck et al. [43].

Climate Types	A (Tropical)	B (Arid)	C (Temperate)
	D (Cold)	E (Polar)
Land Cover Types	Evergreen Needleleaf Forest (ENF)	Evergreen Broadleaf Forest (EBF)	Deciduous Needleleaf Forest (DNF)
	Deciduous Broadleaf Forest (DBF)	Mixed Forest (MIF)	Open Shrublands (OSH)
	Woody Savannas, Savannas (WSA)	Grasslands (GRA)	Croplands (CRO)
	Cropland and Natural Vegetation Mosaic (CNV)	Snow and Ice (SNI)	Barren or Sparsely Vegetated (BSV)

2.4. Statistical Analysis Strategy

For each pixel, WY was estimated as follows:

$$\mathrm{WY}=\mathrm{P}-\mathrm{ET},$$

(1)

where WY, P, and ET are the annual water yield, annual precipitation, and annual ET in mm/a, respectively. We used the MK trend test to reduce the effect of outliers [44]. This is commonly used as a test of significance for the Theil–Sen estimator, which is a nonparametric method to estimate a monotonic trend from a time series. Temporal autocorrelation was accounted for by AR correction. For a time series, each value is not only related to the previous values, but also the resulting disturbance. Auto-regressive (AR) model is a statistical correction approach to time series. Not accounting for temporal autocorrelation could lead to increased false-positive rates. The procedure outlined by von Storch [45] was followed. The lag operator (lag-1) is an operator that converts the previous value of a time series into the current value. Von Storch [45] pointed out that when performing the MK test, the “pre-whitened” time series was much less troubled by serial correlation. That is, the lag-1 autocorrelation is estimated first, and then the original time series is replaced. For each pixel, the temporal autocorrelation was calculated at lag-1 as follows:

$$\widehat{\mathrm{r}}=\frac{\mathrm{n}·{{\displaystyle \sum }}_1^{\mathrm{n}-1}\left({\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{x}}\right)\left({\mathrm{x}}_{\mathrm{i}+1}-\overline{\mathrm{x}}\right)}{\left(\mathrm{n}-1\right)·{{\displaystyle \sum }}_1^{\mathrm{n}}{\left({\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{x}}\right)}^2},$$

(2)

where n was the number of years, ${\mathrm{x}}_{\mathrm{i}}$ was the data for year i, and $\overline{\mathrm{x}}$ was the multi-year average. After that, the original time-series, ${\mathrm{x}}_{\mathrm{i}}$, was replaced with the series, ${\mathrm{y}}_{\mathrm{i}}$ as follows:

$${\mathrm{y}}_{\mathrm{i}}={\mathrm{x}}_{\mathrm{i}}-{\widehat{\mathrm{r}}\mathrm{x}}_{\mathrm{i}-1}$$

(3)

The MK trend test was then performed on this new time series.

To correct for multiple hypothesis testing, a clustering-based permutation method was applied (Figure 1). The permutation test established a threshold for overall significance based on the number of contiguous (first-order queen neighbors) significant grid cells. The motivation for this test was that the clusters formed by false positives were smaller than those formed when a true increase or decrease in WY existed. A total of 3000 permutations were performed to determine the threshold for statistical significance. At each permutation, all grid cells were permuted jointly, that is, using the same set of permuted time indices. This preserves the spatial correlation of the data. The MK trend test was performed on the permuted data, and the largest cluster of significant grid cells was recorded. The 2700th of the 3000 largest recorded clusters (i.e., 90th percentile) was the threshold for overall significance. In the original data, the number of grid cells for each significant region was counted, and if it exceeded the threshold established by the permutation method, it was declared that this region had a confirmed change. This method has been proven to control the probability of false positives in the results and has higher statistical power than comparable methods [31].

For the grid cells with a confirmed significant change in WY, the response patterns of WY to precipitation and ET were identified according to the relationships of their Sen’s slopes. Theil–Sen’s slope estimation is usually used to calculate trend values with MK test. That is, Sen’s slope is first calculated, and then the MK test is used to determine the trend significance. This method is computationally efficient, insensitive to measurement errors and outlier data, and is often used in trend analysis of long time series data [28,46]. Here, six response patterns were defined as follows: S_WY+, S_P+, and S_ET+ (Pattern Ⅰ); S_WY+, S_P+, and S_ET− (Pattern Ⅱ); S_WY+, S_P−, and S_ET− (Pattern Ⅲ); S_WY−, S_P−, and S_ET− (Pattern Ⅳ); S_WY−, S_P−, and S_ET+ (Pattern Ⅴ); S_WY−, S_P+, and S_ET+ (Pattern Ⅵ), where S_WY+, S_P, and S_ET stand for the Sen’s slope of WY, precipitation, and ET, respectively.

Based on these six patterns, the dominant component leading to the variation in WY was explored. For precipitation and ET varying in the same direction, that with a large Sen’s slope was taken as the dominant factor. When the two were different, the one with a significant sign was considered the dominant factor. If both were significant, precipitation and ET were considered to act together in WY, whereas if neither were significant, they were set to be determined.

3. Results

3.1. Annual WY at the Global Scale

The global 0.5° × 0.5° WY multiyear average (Figure 2a) from 1981 to 2018 was greater than 0 mm in 86.16% of the terrestrial region. In addition, 7.15% of the regions had greater than 1000 mm at 10° S–30° N in eastern Asia, Greenland, and the Amazon basin. The maximum value of 5464.8 mm was located at (77.25°W, 4.25°N). In contrast, 13.84% of the regions had an average WY of less than 0 mm. Among them, 4.82% were less than −500 mm, located on the California Peninsula, the coast of Peru, and the coast of the Arabian Peninsula. The minimum value of −1796.60 mm was located at (115.75°E, 20.25°S). The spatial average of global WY over many years was 299.15 mm/a. The anomalies of 41.59% of the regions with a multi-year average WY greater than 0 mm were greater than 0 mm, while those of all regions with a multi-year average WY of less than 0 mm were less than 0 mm (Figure 2b). There were obvious variations in eastern Asia, Europe, the low-and mid-latitudes of Africa, the mid-latitudes of South America, eastern North America, and mid-to high-latitudes.

We obtained the annual trends in global land precipitation, ET, and WY from 1981 to 2018 (Figure 3). The results showed that the Sen’s slopes of annual precipitation, ET, and WY were 0.4401, 0.2634, and −0.1188 mm/a, with p-values of 0.0071, 0.00001, and 0.5462 (obtained by the MK test), respectively. The average annual precipitation, ET, and WY were 645.38, 220.50, and 229.15 mm/a. The maximum value of annual precipitation was 665.15 mm/a in 2011, and the minimum value was 616.68 mm/a in 1987. The maximum values occurred in 1981, 1989, and 1996, while the minimum values occurred in 1987, 1991, and 2001. The maximum value of annual ET was 230.52 mm/a in 2016, and the minimum value was 211.33 mm/a in 1992. The maximum annual WY was 318.32 mm/a in 1996, and the minimum value was 275.10 mm/a in 1987. The maximum difference between annual precipitation and ET was 441.35 mm in 1999, and the minimum difference was 401.86 mm in 1987.

3.2. Inter-Annual Variability of WY

Considering the spatial differences, MK significance tests for each pixel were performed, and 8801 grids with significant changes in global terrestrial WY were obtained (Figure 4a). Not adjusting for multiple hypothesis testing could overestimate the global WY change. Hence, after applying the correction, a confirmed change in 7872 grids was detected (Figure 4b). Areas where there were grids with trend signals that were not strong enough to be detected by the MK method but strong enough to be detected when performing no correction were categorized as uncertain changes. A total of 10.6% of the origin significant changing areas, which are called uncertain changing areas, were screened out after the correction. The mean value of the global terrestrial WY Sen’s slope before the correction was 0.4855, while after the correction, the mean value of the confirmed area was 0.5266 and that of the uncertain area was 0.2496, both showing an overall increasing trend.

Before applying the correction, 3756 grids (5.72%) showed a significant increasing trend, and 5045 grids (7.68%) showed a significant decreasing trend. After correction, 3361 grids (5.12%) of the global terrestrial WY in the confirmed regions showed an increasing trend, and 4511 grids (6.87%) showed a decreasing trend. The confirmed increasing trend of terrestrial WY in each continent was 48.2% in Africa and 31.3% in Asia, whereas the decreasing trend was 36.8% in Asia, 19.4% in North America, 18.6% in Africa, 12.8% in South America, and 11.5% in Australia. In the uncertain areas of terrestrial WY, 395 grids (0.60%) showed an increasing trend, and 534 grids (0.81%) showed a decreasing trend. Among the increasing trend areas, 30.6% were in North America, and 27.2% were in Asia. Decreasing trends were found in Asia (27.0%), North America (25.3%), and Africa (24.8%) (Figure 5).

Based on the climate-plant types, the global terrestrial (Figure 6) WY Sen’s slope was positive in EBF, WSA, and CNV in the tropical (A) climate zone, with a significant decrease in WY in CRO and both an increase and a decrease proportionately in GRA. Moreover, WY decreased in DBF, OSH, CRO, and BSV but increased in CNV, and both increased and decreased proportionately in WSA and GRA in the arid (B) climate zone. In the temperate (C) climate zone, WY decreased in ENF, EBF, MIF, OSH, GRA, and CRO, whereas it increased in WSA and CNV, and WY both increased and decreased proportionately in DBF. In the cold (D) climate zone, ENF, OSH, and GRA showed an increased and decreased WY in proportion to each other, while DNF, DBF, MIF, WSA, CRO, and CNV showed an increased WY. In the polar (E) climate zone, WY increased and decreased in OSH and GRA in proportion to each other, while it showed an overall increase in SNI and a decrease in BSV.

3.3. WY in Response to Variabilities of Its Components

Based on the global terrestrial WY response to precipitation and the ET model (Figure 7), Pattern I (S_WY+, S_P+, and S_ET+) was concentrated in Africa and Russia, and Pattern II (S_WY+, S_P+, and S_ET−) in Africa. Further, Patterns I (59.3%) and II (40.4%) were widely distributed in the areas with an increasing WY. Pattern IV (S_WY−, S_P−, and S_ET−) was scattered around the globe, while Pattern V (S_WY−, S_P−, and S_ET+) was more clearly distributed globally, with some parts of North America and southern South America. Additionally, Pattern VI (S_WY−, S_P+, and S_ET+) was mainly concentrated in the eastern part of Africa.

Based on the climate-plant types, the WY response pattern to the variabilities of its components was obtained (Figure 8). Areas in A-CRO, B-CRO, and C-CRO with Pattern V accounted for the majority, whereas the D-CRO was associated with Pattern I. In the MIF of different climate zones, different patterns occupied the following different proportions: C-MIF with Pattern V, D-MIF with Pattern I, and E-MIF with Pattern V. The difference was that the CNV and OSH in different regions showed the same pattern.

3.4. Dominating Factors for the Changes in WY

The spatial distribution of the dominant factors of global WY change (Figure 9) showed that the area of the region where WY changes were dominated by precipitation accounted for 31.5% and was distributed in central Africa, North America, and Siberia. The areas dominated by a combination of precipitation and ET accounted for 61.2% and were distributed worldwide. The regions where WY changes were dominated by ET were found on the western side of Africa. The regions with a significant increasing trend in WY were dominated by precipitation (59.2%) and the combined effect (40.4%), while the regions with a significant decreasing trend were mainly dominated by the combined effect (78.4%).

Based on the climate-plant types, the main factors of WY variation were mapped (Figure 10). The WY changes in regions with A-CRO, B-CRO, and C-CRO were dominated by the combined effect, whereas D-CRO was dominated by precipitation. Similarly, C-MIF and E-MIF were dominated by the combined effect, whereas D-MIF was dominated by precipitation. The areas with A-CNV and C-CNV were dominated by precipitation and combined effects, whereas the areas with B-CNV and D-CNV were dominated by precipitation. The OSH in all climate regions was primarily a combined effect.

4. Discussion

Improving the understanding of changes in WY and the dominant drivers is essential for assessing the impact of global climate change on the water cycle. Previous studies have suggested that WY changes due to precipitation and ET in local areas. However, spatiotemporal changes in WY in global terrines and their response patterns to precipitation and ET have not yet been considered. In this study, a new step, correcting for multiple hypothesis testing, was introduced. This study showed that the increasing and decreasing trends in WY have changed around the world. Furthermore, the contribution of the WY change was different in different regions.

The high values of the multi-year average WY were mostly concentrated in coastal areas, while the low values were concentrated inland (Figure 2a). As for the trend of “Dry gets drier, wet gets wetter”, we found (Figure 2b) that WY in all dry regions where multi-year average WY was less than 0 mm/a was decreasing. It was a negative indicator of local water resources [47]. However, for areas that are already rich in WY, it is still necessary to pay attention to the climate change and land cover change to avoid the impact of the decreasing WY somewhere. Some researchers [48] also pointed out that the impact of climate change on global WY was reflected in wet and dry seasons.

The application of the correction based on clustering corrected the overestimation of the degree of global WY change (Figure 4). After applying the correction, we found that approximately 10% of the global areas that demonstrated a significant change in uncertainty were eliminated. The extent of overestimation can be better understood by examining climate zones and land-use types. The results showed that the uncertain areas were mainly concentrated in non-woody plants, such as open shrublands, closed shrublands, and savannas. Moreover, the number of uncertain grids was relatively high for arid desert, arid steppe, and polar tundra climates. In addition, harsher climates and sparser vegetation were observed. WY was mostly increasing in the tropical zone (15°S–15°N), while an overall decreasing trend was observed near 30°N and 30°S (the global desert strip), respectively. We also found that the area in the arid climate type appeared to have more decreasing WY trends. In addition, it is worth noting that the changes were not symmetrical in the northern and southern hemispheres. At high latitudes in the Northern Hemisphere, both increasing and decreasing trends were observed (Figure 11). These might relate to the global wind belts.

Global terrestrial WY during 1981–2018 showed a fluctuating downward trend (Figure 3). Based on the latitudinal variation map of global terrestrial WY changes (Figure 4), WY at low latitudes was dominated by an increasing trend, that at mid-latitudes in the northern and southern hemispheres was dominated by a decreasing trend, and that at high latitudes in the northern hemisphere had a proportion of increasing and decreasing trends. We selected some typical regions for the discussion (Figure 12). All regions showed a consistent WY increasing or decreasing trend, except central Africa (Figure 13). According to our defined changing patterns and dominators, Canada (Figure 12(ⅰ)) showed a decreasing trend in WY, caused by a decrease in precipitation and ET. Here, the multi-year average WY was 95.96 mm/a, and Sen’s slope was −1.65 mm/a. Li et al. [49] also showed that the most obvious decreasing trend of WY in Canada occurred in the northeastern and southern mountains but was attributed to the decrease in precipitation and increase in ET. The influence of glacier/permanent snowmelt from the Rocky Mountains and permafrost degradation could be the potential factors in Northwest Canada because of global warming [50,51]. In contrast, the decrease in WY in West Asia (Figure 12(ⅱ)) was mainly caused by a combination of decreased precipitation and increased ET. There was a quite low multi-year average WY (25.17 mm/a) where Sen’s slope was −1.57 mm/a. The combination of decreasing precipitation and increasing atmospheric demand led to a significant increase in aridity in many subtropical areas [52]. Shirmohammadi et al. [53] suggested that local land-use changes might have a higher impact on WY than climate change. The conversion of rangeland and farmland to orchards would decrease WY. These land-use changes affect WY by a structural change in the basin due to an increase in infiltration and ET rates [54]. An overall drying trend in Southern Hemisphere semiarid regions has been noted (since the 1950s for SE Australia [55]). In our results, the decrease in WY (Sen’s slope was −2.16 mm/a) was due to decreased precipitation over most of southeastern Australia (Figure 12(ⅲ)). The substantially higher temperatures, which enhanced atmospheric evaporative demand, coupled with below-average precipitation, resulted in a strong region-wide drought throughout the entire southeastern Australia [56]. Moreover, it might be influenced by sea surface temperature and its interactions with the atmosphere [57].

The increase in WY (Sen’s slope was 2.8573 mm/a) in Siberia (Figure 12(ⅳ)) was caused by an increase in precipitation, with a multi-year average WY of 323.94 mm/a. Extra moisture could be redistributed to the surface through large-scale circulation, creating precipitation in downwind regions. Lian et al. [58] suggested that spring moisture input from the strongly greened upstream regions (Europe) with westerly winds helped to maintain positive feedback of summer precipitation in Siberia when it became weaker. This notable teleconnection resulted in a localized increased wetness trend. Remote correlation is also applied to the Amazon region (Figure 12(ⅴ)). In our study, the Amazon region was mostly characterized by an increase in WY due to increased precipitation and a combination of increased precipitation and decreased ET. The multi-year average WY was high (684.60 mm/a), and Sen’s slope was 7.87 mm/a. Warming from the tropical Atlantic during the wet season triggered more latitudinal water vapor. Seasonal precipitation was important for the maintenance of tropical rainforest ecosystems. Seager et al. [59] suggested that the changes in WY induced by rising specific humidity, which accompanied atmospheric warming, were offset by the slowdown of the tropical circulation to some extent. The large area of the 10° N east-west strip-trending region of Africa (Figure 12(ⅵ)) showed a positive value of Sen’s slope in WY, with a multi-year average WY of 276.92 mm/a, mostly due to increased precipitation. The increase in WY was partly caused by an increase in precipitation and a decrease in ET. Possibly, the increase in WY was due to greening locally or an increase in precipitation in the basin. Data in the study by Cortés et al. [60] showed a significant trend of greening in the area north of the equator in Central Africa. On one hand, some studies suggested that deforestation could increase WY [21], while greening decreased WY by increasing ET [61,62]. On the other hand, Makarieva et al. [63] argued that increasing ET could drive water vapor transport by changing atmospheric pressure. The amount of water recovered to the ground depended greatly on the watershed area that is, the larger the geographical area, the higher the recovery potential.

This study had some limitations. First, the land cover change affected the land surface water cycle, ET, and precipitation [64]. In our study, the IGBP for 2001, in the middle of the study period, was selected as the classification map of the surface cover to minimize the impact of land cover changes on data monitoring. Balist et al. [65] studied the changes in the WY of the Sirvan River Basin in Iran from 2013 to 2019. They found that forests decreased and built-up areas increased. Simultaneously, the precipitation decreased, and temperatures increased in this area over a long period. Therefore, they concluded that the effect of these factors on the WY is a complex process. Second, because we only focused on monotonic trends for the whole period, segmentation might have occurred (e.g., increasing first, then decreasing), which led to overall nonlinear changes. Testing such segmented trends while considering multiple hypothesis testing remains to be investigated. Third, some regions were too small to be detected in our study. Therefore, analysis of these smaller areas might reveal previously undetected trends. Fourth, we found that significant clustering surrounds some insignificant grids. Noise in the trend might lead to insignificance because there was no wall-to-wall ground-truth validation. Nevertheless, our study explicitly showed that global terrestrial WY variability is influenced by precipitation and ET, with different regions showing corresponding trends of increasing or decreasing WY. This study provides a reference for how WY variability can be used for local water management in the future.

5. Conclusions

In summary, this study revisited significant changes in the terrestrial WY using an improved MK test and investigated the response patterns of WY to precipitation and ET. The traditional MK test showed that the terrestrial WY changed significantly over 13.4% of continents, whereas 10.6% of which were removed because of their uncertainties in the permutation test. The confirmed increasing trends of WY were located in Africa, eastern North America, and Siberia, whereas the confirmed decreasing trends were concentrated in Asia and Oceania. The dominant factors for the increase in WY were mainly precipitation, while the combined effect of precipitation and ET accounted for most of the significant decreases in WY.