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

: Water yield (WY) refers to the difference between precipitation and evapotranspiration (ET), which is vital for available terrestrial water. Climate change has led to signiﬁcant 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 signiﬁcant ( p < 0.05) increase of 5.72% of the total valid grid cells, while it showed a signiﬁcant decrease of 7.68% of those. After correction using the calibration method, the signiﬁcantly increasing and decreasing areas are reduced by 10.52% and 10.58% of them, respectively. After the correction, the conﬁrmed 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 beneﬁcial for improving the understanding of WY in response to hydrological variables in the context of climate change. and 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.


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

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].

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: Remote Sens. 2022, 14,2009 3 of 17 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.

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).

Statistical Analysis Strategy
For each pixel, WY was estimated as follows: 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:r where n was the number of years, x i was the data for year i, and x was the multi-year average. After that, the original time-series, x i , was replaced with the series, y i as follows: 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]. 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.

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 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 I); S WY +, S P +, and S ET − (Pattern II); S WY +, S P −, and S ET − (Pattern III); S WY −, S P −, and S ET − (Pattern IV); S WY −, S P −, and S ET + (Pattern V); S WY −, S P +, and S ET + (Pattern VI), 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.

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.

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. 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 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.
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.

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 (SWY+, SP+, and SET+) was concentrated in Africa and Russia, and Pattern II (SWY+, SP+, and SET−) in Africa. Further, Patterns I (59.3%) and II (40.4%) were widely distributed in the areas with an increasing WY. Pattern IV (SWY−, SP−, and SET−) was scattered around the globe, while Pattern V (SWY−, SP−, and SET+) was more clearly distributed globally, with some parts of North America and southern South America. Additionally, Pattern VI (SWY−, SP+, and SET+) 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 propor-

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.     Figure 6) across the grid cells in which WY significantly changed during 1981-2018 after correcting for multiple hypothesis testing (α = 0.05). SWY, SP, and SET indicate the Theil-Sen's slope of WY, P, and ET, respectively. The symbols + and − indicate a positive and negative WY (P or ET) trend during the study period, respectively. Figure 6) across the grid cells in which WY significantly changed during 1981-2018 after correcting for multiple hypothesis testing (α = 0.05). S WY , S P , and S ET indicate the Theil-Sen's slope of WY, P, and ET, respectively. The symbols + and − indicate a positive and negative WY (P or ET) trend during the study period, respectively.

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%).

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.

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.

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(i)) 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(ii)) 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(iii)). 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]. 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]. Figure 11. Significant changes in annual terrestrial water yield according to the number of grid cells in different altitude. Darker shades indicate the confirmed increasing (decreasing) trends that remains after correcting for multiple hypothesis testing (α = 0.05), while lighter shades indicate the uncertain trends, i.e., additional significant changes that are detected when no multiple testing correction took place.
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 Figure 11. Significant changes in annual terrestrial water yield according to the number of grid cells in different altitude. Darker shades indicate the confirmed increasing (decreasing) trends that remains after correcting for multiple hypothesis testing (α = 0.05), while lighter shades indicate the uncertain trends, i.e., additional significant changes that are detected when no multiple testing correction took place.
The increase in WY (Sen's slope was 2.8573 mm/a) in Siberia (Figure 12(iv)) 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(v)). 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(vi)) 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. 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 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. Figure 12. Statistics for the confirmed trends of terrestrial water yield (WY) represented as the number of grid cells by the latitude. Blue and red indicate a significant increase and decrease in WY during 1981-2018 after correcting for multiple hypothesis testing (α = 0.05), respectively. The red boxes indicate the areas of interest, including parts of Canada (ⅰ), Siberia (ⅱ), West Asia (ⅲ), Australia (ⅳ), Central Africa (ⅴ), and the Amazon (ⅵ). 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 Figure 13. Temporal variabilities of ten-year average terrestrial water yield (WY) from 1981 to 2018 in typical regions ( Figure 12). The colored boxes indicate the areas of interest, including parts of Canada (i), West Asia (ii), Australia (iii), Siberia (iv), the Amazon (v), and Central Africa (vi).

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.