Estimating Rainfall Interception of Vegetation Canopy from MODIS Imageries in Southern China

The interception of rainfall by vegetation canopies plays an important role in the hydrologic process of ecosystems. Most estimates of canopy rainfall interception in present studies are mainly through field observations at the plot region. However, it is difficult, yet important, to map the regional rainfall interception by vegetation canopy at a larger scale, especially in the southern rainy areas of China. To obtain a better understanding of the spatiotemporal variation of vegetation canopy rainfall interception with regard to the basin scale in this region, we extended a rainfall interception model by combining the observed rainfall data and moderate resolution imaging spectroradiometer leaf area index (MODIS_LAI) data to quantitatively estimate the vegetation canopy rainfall interception rate (CRIR) at small/medium basin scales in Guangdong Province, which is undergoing large changes in vegetation cover due to rapid urban expansion in the area. The results showed that the CRIR in Guangdong declined continuously during 2004–2012, but increased slightly in 2016, and the spatial variability of CRIR showed a diminishing yearly trend. The CRIR also exhibited a distinctive spatial pattern, with a higher rate to the east and west of the mountainous areas and a lower rate in the central mountainous and coastal areas. This pattern was more closely related to the spatial variation of the LAI than that of rainfall due to frequent extreme rainfall events saturating vegetation leaves. Further analysis demonstrated that forest coverage, instead of background climate, has a certain impact on the canopy rainfall interception, especially the proportion of broad-leaved forests in the basin, but more in-depth study is warranted in the future. In conclusion, the results of this study provide insights into the spatiotemporal variation of canopy rainfall interception at the basin scale of the Guangdong Province, and suggest that forest cover should be increased by adjusting the species composition to increase the proportion of native broad-leaved species based on the local condition within the basin. In addition, these results would be helpful in accurately assessing the impacts of forest ecosystems on regional water cycling, and provide scientific and practical implications for water resources management.


Introduction
Canopy rainfall interception is the amount of rainwater that is intercepted, stored, and subsequently lost by evaporation from the canopy [1,2].It is a significant component of the water cycle in the vegetation ecosystem since it influences many hydrological processes such as infiltration, erosion, soil moisture distribution, runoff dynamics, and flood generation [3,4].Considering that canopy rainfall interception is an important component of the water balance in a terrestrial ecosystem, an accurate estimation of rainfall interception by the vegetation canopy is of great value in the study of hydrological processes [5].This hydrological process is also very important in water resource management and in the context of climate change [2,6].
Since canopy rainfall interception cannot be measured directly, it is usually estimated by determining the difference between gross rainfall and the sum of throughfall and stemflow [7][8][9], which are commonly assessed through the study of individual leaves or whole plants and extrapolated to entire plant communities.This means that an average value is therefore generally obtained by using a large number of collectors and a large collection area [10,11].Ultimately, these methods are time-consuming and expensive to implement.Although these field observations provide high precision data on a small scale, the extrapolation of rainfall interception data to a larger area is difficult to achieve, and the spatiotemporal patterns of interception are also difficult to map quantitatively.Many scientists have attempted to find a cost-effective approach in documenting changes over large geographic regions in both time and space [2,[12][13][14].Modeling has become an effective tool in understanding the spatial distribution of canopy interception at a larger scale [6,15,16].Many models of the rainfall interception process have been developed [17][18][19].The best-known models for canopy rainfall interception are the Rutter model [20,21] and Gash model [22].One of the most important parameters in the Gash and Rutter models is canopy storage capacity (S max ) [22][23][24].This parameter indicates that exceeding the threshold storm size most likely results in a significant overflow of water through the canopy and thus downward toward the soil surface [5,25].Therefore, canopy storage capacity is a key factor in controlling rainfall interception by vegetation canopy [7,26,27].
However, estimating an average value of S max in large or medium-sized areas is difficult due to the high spatial variability of vegetation.Many studies have demonstrated that canopy rainfall storage capacity is generally related to the leaf area index (LAI), which can be derived from remote-sensing images [28,29].Many publications have shown that remote sensing is widely recognized as an effective tool for spatially and temporally monitoring the changing patterns of vegetation, and it has significant advantages over in situ measurements [6,30].The ability to rapidly assess the LAI using vegetation indices from remote-sensing images provides a means of determining the phenological characteristics of cover vegetation over a wide geographic area [31][32][33].
Guangdong Province is situated in a subtropical zone in southern China, which has been confronted by a large change in vegetation cover due to the rapidly developing economy and subsequent urban expansion [34,35].In addition, many studies have shown that, in the context of climate change, the number of extreme rainfall days and the intensity of storms have increased, and abnormal rainfall and flood emergencies have constantly occurred in Guangdong [36,37].For these reasons, in addition to its unique soil characteristics with thin soil layers, disproportionate soil texture, and its specific topography and landform characteristics with mountain and hilly areas occupying approximately 60 % of the total area [38,39], Guangdong is particularly susceptible to soil erosion, especially that caused by rainfall erosive processes [40].Given that rainfall interception by vegetation canopy, which accounts for a large proportion of atmospheric precipitation (totals were between 10% and 50% of annual gross rainfall) [41][42][43], plays an important role in ecosystem hydrological processes including infiltration, erosion, and even runoff [44], it is all the more important to quantitatively simulate and analyze vegetation canopy rainfall interception as well as its spatiotemporal variation characteristics in Guangdong Province.However, studies have been focused more frequently on the characteristics of rainfall interception at the plot scale and less often at the basin scale because impact factors such as precipitation and vegetation characteristics are difficult to model spatially.
In this study, a site-based atmospheric rainfall interception of the vegetation canopy model was applied to moderate resolution imaging spectroradiometer leaf area index (MODIS_LAI) remote sensing data coupled with observed precipitation data from weather stations.We quantitatively estimated the rainfall interception capacity of the vegetation canopy in Guangdong at the small/medium basin scale from 2004 to 2016.This is a novel application of remote sensing technology in a combination of ecological models.The results could be very helpful in accurately assessing the impacts of forest ecosystem on regional water cycling and provide scientific and practical implications for water resource management.

Study Area
Guangdong Province is located in southern China (20 • 14 -25 •  The soil distribution in Guangdong is characterized by an obvious zonality.Hapic Acrisol, Haplic Ferralsol, and Rhodic Ferralsol are distributed from north to south in this region, which accounts for 37.96%, 24.8%, and 5.15% of soil composition, respectively.The tropical and subtropical broad-leaved evergreen forests are the major climax forest communities.The sclerophyllous evergreen savanna shrubs and grasslands are the secondary climax forest communities [46].However, due to poor soil fertility and interference by humans, the original forest community structure has been severely damaged.To reduce the environmental problems of the areas, forestry protection policies and afforestation projects have been implemented, and the amount of forestland has since increased [46,47].However, the number of natural forests has decreased while the number of planted forests has increased. The main landform type in Guangdong is mountainous, and between the mountains there are different sized valleys and basins.Guangdong is divided into 1367 relatively independent and closed small basins, and the drainage area of each small basin is lower than 50 km 2 (Figure 1).In addition, seven major basins, with drainage areas larger than 3000 km 2 , were selected to be the main research objects of this study including the Hanjiang (HJ) basin, Dongjiang (DJ) basin, Beijiang (BJ) basin, Xijiang (XJ) basin, the Pearl River (PR) basin, Moyangjiang (MYJ) basin, and Jianjiang (JJ) basin (Figure 1).The regions beyond these seven major basins with incomplete drainage areas were classified as 'other areas' in this study.impacts of forest ecosystem on regional water cycling and provide scientific and practical implications for water resource management.

Study Area
Guangdong Province is located in southern China (20°14′-25°31′N, 109°40′-117°20′E) and occupies an area of 179,752 km 2 .The climate of Guangdong is tropical and subtropical monsoon and the 12 months can be divided into humid (April to September) and dry (October to March) seasons.The annual mean temperature is 22 °C, and the annual precipitation amount is about 1770 mm, of which 1384 mm (78%) and 386 mm (22%) fall in the wet and dry seasons, respectively [45].
The soil distribution in Guangdong is characterized by an obvious zonality.Hapic Acrisol, Haplic Ferralsol, and Rhodic Ferralsol are distributed from north to south in this region, which accounts for 37.96%, 24.8%, and 5.15% of soil composition, respectively.The tropical and subtropical broad-leaved evergreen forests are the major climax forest communities.The sclerophyllous evergreen savanna shrubs and grasslands are the secondary climax forest communities [46].However, due to poor soil fertility and interference by humans, the original forest community structure has been severely damaged.To reduce the environmental problems of the areas, forestry protection policies and afforestation projects have been implemented, and the amount of forestland has since increased [46,47].However, the number of natural forests has decreased while the number of planted forests has increased.
The main landform type in Guangdong is mountainous, and between the mountains there are different sized valleys and basins.Guangdong is divided into 1367 relatively independent and closed small basins, and the drainage area of each small basin is lower than 50 km 2 (Figure 1).In addition, seven major basins, with drainage areas larger than 3000 km 2 , were selected to be the main research objects of this study including the Hanjiang (HJ) basin, Dongjiang (DJ) basin, Beijiang (BJ) basin, Xijiang (XJ) basin, the Pearl River (PR) basin, Moyangjiang (MYJ) basin, and Jianjiang (JJ) basin (Figure 1).The regions beyond these seven major basins with incomplete drainage areas were classified as 'other areas' in this study.

Rainfall Data
Based on ArcGIS software, the gridded daily time-series images of rainfall data with a resolution of 500 m were produced by applying the interpolation method to remain consistent with the moderate resolution imaging spectroradiometer leaf area index (MODIS_LAI) gridded data of 500 m resolution.Three interpolation methods (e.g., inverse distance weighting (IDW), Kriging, and Spline) were fully evaluated before we chose a suitable method to use in this study.A point-by-point assessment was conducted and involved a comparison of the interpolated daily rainfall data in the grid cell of 500 m resolution from June to August 2015 by reference to the 238 situ-observed rainfall data in Wu et al. (2019) [48].The rainfall data were ultimately interpolated to grid cells by using the IDW method due to its higher accuracy in this study (Figure S1).These grid cell data were interpolated with daily precipitation data from January 2004 to December 2016 collected from 87 in situ observation stations belonging to the Guangdong Provincial Meteorological Bureau.

Leaf Area Index (LAI) Data
Due to the low resolution of the early remote sensing data and the uncertainty of the products' data quality from the different sources, it is difficult to ensure the comparability of the products.In this study, the MODIS terrestrial standard product MOD13A1 (16-day composited vegetation index product with a spatial resolution of 500 m) was acquired between February 2004 and December 2016 in hierarchical data format (HDF) from the geospatial data cloud [49], which is maintained by the Scientific Data Center of the Computer Network Information Center, Chinese Academy of Sciences.
Compared with the National Oceanic and Atmospheric Administration (NOAA) satellite and other Landsat data, the rate of space separation has significantly improved.The spatial resolution has increased by one order of magnitude from the kilometer level of NOAA to the 100-meter level of MODIS.Second, the time resolution is more advantageous.The MODIS satellite passes through four times a day (while the NOAA satellite passes through only two times a day) and has a greater real-time monitoring capacity for a variety of sudden, rapidly changing natural disasters.Third, the spectral resolution has been greatly improved.The MODIS satellite has 36 wavebands (while the NOAA satellite only has five wavebands), and its multi-channel observation enhances the ability to observe the Earth's complex systems and surface types.
The MODIS_LAI data were spliced and clipped by ENVI software, and the LAI data of Guangdong were extracted.Next, the data were averaged by ArcGIS software to obtain the monthly and annual LAI gridded data with a resolution of 500 m for 2004 to 2016.In order to ensure the accuracy of the results, this study only selected the high-quality MODIS_LAI data for analysis.
In order to more clearly reveal the spatial and temporal variation patterns in Guangdong during the period of 2004 to 2016, we only present the results and spatial quantization maps for the years 2004, 2008, 2012, and 2016 in the text.The integrated spatial quantization maps from 2004 to 2016 are presented in the Supplementary Materials.

Canopy Rainfall Interception Model
Rainfall interception by vegetation canopy is simulated by calculating a maximum storage capacity (mm, S max ), which is filled during rainfall.The maximum interception storage capacity is estimated using an equation developed by De Roo et al. [50]: where LAI is the leaf area index.
Cumulative interception during each rainfall event is simulated using an equation developed by Aston [51], which is modified from Merriam [52]: where S v is the canopy rainfall interception (mm); P cum is the cumulative rainfall (mm) during each rainfall event; and η is the correction factor, equal to (0.046 × LAI).The theoretical assumption of the equation is when P cum = 0, S v = 0, and P cum → ∞ , S v = Maximum storage capacity (S max ).

Canopy Rainfall Interception Rate (CRIR)
The CRIR of vegetation is calculated by the following equation: where a is the CRIR (%).
To verify the precision and validity of the interception model, the canopy rainfall interception rate of 24 typical and relative heavy precipitation events from our previous study [48] were chosen as the observed data to validate the corresponding modeled rainfall interception rate.Through the verification of error analysis, the accuracy and reliability of the model was feasible in this study (Figure S2).

Propensity Score Analysis
To clarify the temporal variation characteristics of each basin, a unary linear regression model was used to calculate the change in slope from 2004 to 2016, namely the propensity score (SLOPE), which was used to analyze the linear trend of CRIR for each basin [53]: where n is the total number of years (13); x i is the i year (2004 is the first year); and A i represents the corresponding vegetation CRIR for the i year.

Categorization of Canopy Rainfall Interception
To reflect the spatial distribution characteristics of the vegetation canopy rainfall interception in Guangdong, based on the canopy rainfall interception of all small basins in 2004, the subregions of Guangdong were divided into three areas with different levels of interception (high, moderate, and low), as shown in Table 1.

Spatiotemporal Variation Pattern of Rainfall
The annual rainfall of Guangdong in 2004, 2008, 2012, and 2016 was 1334.18 mm, 1531.58 mm, 1875.43 mm, and 2327.89mm, respectively.The overall rainfall in Guangdong has increased due to the impacts of climate change [44].The most noticeable increases in rainfall occur along the coast and in the eastern part of Guangdong (Figure 2, Figure S3), which are mainly caused by the increase in the number of extreme weather events such as subtropical cyclones.In general, the amount of rainfall in Guangdong mainly presents a spatial pattern of decreasing gradually from south to north, and this pattern is particularly prominent in wetter years such as 2012 and 2016 (Figure 2).
For the seven major basins (Figure 3 and Figure S4), the rainfall trend for each basin was the same as the overall increasing rainfall trend in Guangdong.Specifically, the annual growth rate of rainfall in the XJ and PR basins is particularly noticeable, with a growth range between 0 and 1000 mm.The rainfall in the BJ, DJ, and HJ basins slightly increased from 2004 to 2008, with a growth range between 0 and 90 mm, but then rapidly increased from 2008 to 2016, with a growth range between 0 and 1200 mm.In the MYJ and JJ basins, the annual rainfall rapidly increased from 2004 to 2008, gradually increased from 2008 to 2012, and then slightly decreased from 2012 to 2016.In addition, the annual mean rainfall displayed obvious zonality, which was significantly greater in coastal regions than in inland regions.Areas with relatively little rainfall were mostly found in the basin with high surroundings and a low middle point (e.g., the BJ basin).However, areas with relatively high rainfall were generally distributed across flat coastal areas or in low-lying valleys (e.g., the MYJ basin).For the seven major basins (Figures 3 and S4), the rainfall trend for each basin was the same as the overall increasing rainfall trend in Guangdong.Specifically, the annual growth rate of rainfall in the XJ and PR basins is particularly noticeable, with a growth range between 0 and 1000 mm.The rainfall in the BJ, DJ, and HJ basins slightly increased from 2004 to 2008, with a growth range between 0 and 90 mm, but then rapidly increased from 2008 to 2016, with a growth range between 0 and 1200 mm.In the MYJ and JJ basins, the annual rainfall rapidly increased from 2004 to 2008, gradually increased from 2008 to 2012, and then slightly decreased from 2012 to 2016.In addition, the annual mean rainfall displayed obvious zonality, which was significantly greater in coastal regions than in inland regions.Areas with relatively little rainfall were mostly found in the basin with high surroundings and a low middle point (e.g., the BJ basin).However, areas with relatively high rainfall were generally distributed across flat coastal areas or in low-lying valleys (e.g., the MYJ basin).

Spatiotemporal Variation Patterns of LAI
From Figures 4 and S5, we found that the spatial heterogeneity of the land surface vegetation is diverse in Guangdong.The annual mean LAI was relatively low in the central Pearl River Delta region, the southwest Leizhou Peninsula, and the east Chaoshan Plain in Guangdong.This is mainly due to the rapid economic development and exploitation of natural resources in the Pearl River Delta region, and the perennial typhoon in the Chaoshan Plain.In addition, the Leizhou Peninsula's terrain gradually becomes lower from north to south, and the streamflow in its short, shallow river quickly drains into the sea.Therefore, it is difficult to build a reservoir there to store water, which makes the peninsula prone to drought.In addition, since the peninsula is located in a subtropical monsoon region, typhoons also affect this region from May to October every year.The above reasons lead to a poor growth environment for vegetation and results in a lower LAI value for the region.The annual mean LAIs in the east, west, and north mountainous areas were higher than that in other areas of Guangdong (Figure 4).The vegetation coverage and quality were correspondingly higher in these regions, which is mainly due to better vegetation sites and growing conditions.The ecological background such as topography and geomorphology as well as the degree of economic development are important factors affecting the spatial distribution patterns of vegetation in Guangdong.The annual mean LAIs of this area in 2004, 2008, 2012, and 2016 were 1.50, 1.57, 1.59, and 1.80, respectively, which indicates that the overall vegetation coverage in Guangdong has increased (Figure S5).
According to the analysis on the seven major river basins (Figure 5), the annual mean increase in the LAI in the BJ, XJ, and MYJ basins were the largest, increasing from 1.61, 1.51, and 1.69 in 2004 to 2.09, 1.97, and 2.14 in 2016, respectively.In addition, the annual mean LAI in the PR and DJ basins increased slightly from 0.98 and 1.71 in 2004 to 1.27 and 1.91 in 2016, respectively.However, the annual increase of the LAI in the HJ basin slightly increased from 2004 to 2012, but decreased after 2012 (Figure S6).

Spatiotemporal Variation Patterns of LAI
From Figure 4 and Figure S5, we found that the spatial heterogeneity of the land surface vegetation is diverse in Guangdong.The annual mean LAI was relatively low in the central Pearl River Delta region, the southwest Leizhou Peninsula, and the east Chaoshan Plain in Guangdong.This is mainly due to the rapid economic development and exploitation of natural resources in the Pearl River Delta region, and the perennial typhoon in the Chaoshan Plain.In addition, the Leizhou Peninsula's terrain gradually becomes lower from north to south, and the streamflow in its short, shallow river quickly drains into the sea.Therefore, it is difficult to build a reservoir there to store water, which makes the peninsula prone to drought.In addition, since the peninsula is located in a subtropical monsoon region, typhoons also affect this region from May to October every year.The above reasons lead to a poor growth environment for vegetation and results in a lower LAI value for the region.The annual mean LAIs in the east, west, and north mountainous areas were higher than that in other areas of Guangdong (Figure 4).The vegetation coverage and quality were correspondingly higher in these regions, which is mainly due to better vegetation sites and growing conditions.The ecological background such as topography and geomorphology as well as the degree of economic development are important factors affecting the spatial distribution patterns of vegetation in Guangdong.The annual mean LAIs of this area in 2004, 2008, 2012, and 2016 were 1.50, 1.57, 1.59, and 1.80, respectively, which indicates that the overall vegetation coverage in Guangdong has increased (Figure S5).
According to the analysis on the seven major river basins (Figure 5), the annual mean increase in the LAI in the BJ, XJ, and MYJ basins were the largest, increasing from 1.61, 1.51, and 1.69 in 2004 to 2.09, 1.97, and 2.14 in 2016, respectively.In addition, the annual mean LAI in the PR and DJ basins increased slightly from 0.98 and 1.71 in 2004 to 1.27 and 1.91 in 2016, respectively.However, the annual increase of the LAI in the HJ basin slightly increased from 2004 to 2012, but decreased after 2012 (Figure S6).

Variation Pattern of Annual Canopy Rainfall Interception
Based on the cumulative rainfall canopy interception in 2004, the small basins in Guangdong were divided into three canopy rainfall interception grades (high, moderate, and low interception area, Figure 6) according to the standard categorization of rainfall canopy interception described in Table 1.
In Table 2, we can see that the annual rainfall canopy interception for the different regions gradually increased over the years.The annual rainfall canopy interception in the moderate interception area was 131.83 mm in 2016, which exceeded its threshold of 90 mm.The same trend also occurred in the low interception area, with an annual canopy rainfall interception of 70.26 mm in 2016, which also exceeded its threshold of 60 mm.This means that when the annual rainfall and vegetation cover increased, the number of high-magnitude interception areas (e.g., moderate and high interception areas) also increased.Based on the cumulative rainfall canopy interception in 2004, the small basins in Guangdong were divided into three canopy rainfall interception grades (high, moderate, and low interception area, Figure 6) according to the standard categorization of rainfall canopy interception described in Table 1.
In Table 2, we can see that the annual rainfall canopy interception for the different regions gradually increased over the years.The annual rainfall canopy interception in the moderate interception area was 131.83 mm in 2016, which exceeded its threshold of 90 mm.The same trend also occurred in the low interception area, with an annual canopy rainfall interception of 70.26 mm in 2016, which also exceeded its threshold of 60 mm.This means that when the annual rainfall and vegetation cover increased, the number of high-magnitude interception areas (e.g., moderate and high interception areas) also increased.The annual mean CRIRs from 2004 to 2016 all showed the same trend where the rates were lowest in March, and began to rise in April, reached their highest values in October, and began to decrease in November until March the following year (Figure 7).In particular, the fluctuations of monthly mean CRIR were depressed through the years.The difference in monthly mean rate (the difference between the highest and the lowest monthly mean CRIRs within a year) decreased from 6.26% in 2004 to 4.41% in 2016.Annual statistical results have shown that the annual mean CRIR had a downward trend in 2004, 2008, and 2012, but increased slightly in 2016 with annual mean values of 4.80%, 4.53%, 4.05%, and 4.13%, respectively (Figure 8).Considering Equations ( 1)-( 3) and the annual

Variation Patterns of Annual CRIR
The annual mean CRIRs from 2004 to 2016 all showed the same trend where the rates were lowest in March, and began to rise in April, reached their highest values in October, and began to decrease in November until March the following year (Figure 7).In particular, the fluctuations of monthly mean CRIR were depressed through the years.The difference in monthly mean rate (the difference between the highest and the lowest monthly mean CRIRs within a year) decreased from 6.26% in 2004 to 4.41% in 2016.Annual statistical results have shown that the annual mean CRIR had a downward trend in 2004, 2008, and 2012, but increased slightly in 2016 with annual mean values of 4.80%, 4.53%, 4.05%, and 4.13%, respectively (Figure 8).Considering Equations ( 1)-( 3) and the annual rainfall and LAI trends in Guangdong, there are two possible reasons for these variations.On one hand, the increase in rainfall will increase the canopy interception to a certain extent, according to Equation (1), as the vegetation leaves have a saturation value for water absorption.Generally speaking, rainfall canopy interception will increase as the amount of rainfall increases, but once the rainfall reaches a certain point, canopy rainfall interception gradually reaches saturation and plateaus.On the other hand, the increase in rainfall can also inhibit canopy interception according to Equation (3).The CRIR decreases if this inhibiting effect is greater than the driving effect of the vegetation LAI increase on the CRIR.rainfall and LAI trends in Guangdong, there are two possible reasons for these variations.On one hand, the increase in rainfall will increase the canopy interception to a certain extent, according to Equation (1), as the vegetation leaves have a saturation value for water absorption.Generally speaking, rainfall canopy interception will increase as the amount of rainfall increases, but once the rainfall reaches a certain point, canopy rainfall interception gradually reaches saturation and plateaus.On the other hand, the increase in rainfall can also inhibit canopy interception according to Equation (3).The CRIR decreases if this inhibiting effect is greater than the driving effect of the vegetation LAI increase on the CRIR.When considering the spatial characteristic (Figures 8 and S7), we found that the CRIR showed a significant spatial difference pattern that was high in the mountainous areas of the eastern and western flanks, but low in the central mountainous areas and the coastal areas in Guangdong in 2004.There were 149 basins with a CRIR greater than 5%, and 1227 basins with a CRIR lower than 5%, accounting for 40.75% and 59.25% of Guangdong's total area, respectively.The maximum and minimum annual mean CRIRs were 11.67% and 0.33%, respectively, in 2004.However, the spatial difference gradually decreased over time.The number of basins with a CRIR greater than 5% decreased to 79 in 2012, accounting for only 10.49% of the total area of Guangdong, and the maximum and minimum annual mean CRIRs were 10.41% and 0.34%, respectively, in 2012.The annual mean CRIR decreased significantly in western and northern Guangdong from 2004 to 2012, and increased slightly from 2012 to 2016.However, the annual mean CRIR in eastern Guangdong, especially in the HJ basin, decreased from 2004 to 2016.Based on the analysis in Sections 3.1 and 3.2, though both the annual mean LAI and annual rainfall gradually increased, the spatial heterogeneity of the annual mean LAI increased (Figure 4) while the spatial heterogeneity of the annual rainfall decreased (Figure 2).Combining the saturation threshold of LAI-Smax (Equation ( 1)) and the inhibitory effect of rainfall on CRIR (Equation ( 3)), the spatial homogeneity of annual rainfall further magnified the spatial pattern of the CRIR, while the spatial heterogeneity of the annual mean LAI narrowed the spatial pattern of the CRIR, leading to a decrease in the spatial variability of CRIR during the 2004-2016 period.From the SLOPE index of CRIR in Guangdong from 2004 to 2016 (Figure 9), we found that the CRIR could potentially increase in the Leizhou Peninsula and west coastal areas of Guangdong.When considering the spatial characteristic (Figure 8 and Figure S7), we found that the CRIR showed a significant spatial difference pattern that was high in the mountainous areas of the eastern and western flanks, but low in the central mountainous areas and the coastal areas in Guangdong in 2004.There were 149 basins with a CRIR greater than 5%, and 1227 basins with a CRIR lower than 5%, accounting for 40.75% and 59.25% of Guangdong's total area, respectively.The maximum and minimum annual mean CRIRs were 11.67% and 0.33%, respectively, in 2004.However, the spatial difference gradually decreased over time.The number of basins with a CRIR greater than 5% decreased to 79 in 2012, accounting for only 10.49% of the total area of Guangdong, and the maximum and minimum annual mean CRIRs were 10.41% and 0.34%, respectively, in 2012.The annual mean CRIR decreased significantly in western and northern Guangdong from 2004 to 2012, and increased slightly from 2012 to 2016.However, the annual mean CRIR in eastern Guangdong, especially in the HJ basin, decreased from 2004 to 2016.Based on the analysis in Sections 3.1 and 3.2, though both the annual mean LAI and annual rainfall gradually increased, the spatial heterogeneity of the annual mean LAI increased (Figure 4) while the spatial heterogeneity of the annual rainfall decreased (Figure 2).Combining the saturation threshold of LAI-S max (Equation ( 1)) and the inhibitory effect of rainfall on CRIR (Equation ( 3)), the spatial homogeneity of annual rainfall further magnified the spatial pattern of the CRIR, while the spatial heterogeneity of the annual mean LAI narrowed the spatial pattern of the CRIR, leading to a decrease in the spatial variability of CRIR during the 2004-2016 period.From the SLOPE index of CRIR in Guangdong from 2004 to 2016 (Figure 9), we found that the CRIR could potentially increase in the Leizhou Peninsula and west coastal areas of Guangdong.From the perspective of the seven major river basins (Figures 10 and S8), the CRIR in the BJ and XJ basins gradually decreased from 2004 to 2012 and slightly increased from 2012 to 2016, with ranges of 0-1.33% and 0-0.74%, respectively.Similarly, the CRIR in the MYJ basin significantly decreased from 2004 to 2008 and then increased slightly from 2008 to 2016.The CRIR in the DJ and HJ basins,   From the perspective of the seven major river basins (Figures 10 and S8), the CRIR in the BJ and XJ basins gradually decreased from 2004 to 2012 and slightly increased from 2012 to 2016, with ranges of 0-1.33% and 0-0.74%, respectively.Similarly, the CRIR in the MYJ basin significantly decreased from 2004 to 2008 and then increased slightly from 2008 to 2016.The CRIR in the DJ and HJ basins, From the perspective of the seven major river basins (Figure 10 and Figure S8), the CRIR in the BJ and XJ basins gradually decreased from 2004 to 2012 and slightly increased from 2012 to 2016, with ranges of 0-1.33% and 0-0.74%, respectively.Similarly, the CRIR in the MYJ basin significantly decreased from 2004 to 2008 and then increased slightly from 2008 to 2016.The CRIR in the DJ and HJ basins, however, decreased from 2004 to 2016 with a range of 0-1.41%.The CRIR in the JJ basin remained essentially invariable during the 2004-2016 period.
however, decreased from 2004 to 2016 with a range of 0-1.41%.The CRIR in the JJ basin remained essentially invariable during the 2004-2016 period.

Effects of Climate and Forest Cover on Variation Patterns of the CRIR
The forest canopy plays an important role in rainfall interception in the forest ecosystem.However, differentiating the effects of climate and forest cover on the spatiotemporal variation patterns of the CRIR is not well understood for this region.To address this issue, we first selected two basins with similar forest coverage (about 70%) in each of the four climatic zones (e.g., middle subtropical zone, southern south subtropical zone, northern south subtropical zone, and north tropical zone, Table 3) in Guangdong to analyze their variation patterns of CRIR in the same climatic zone (Figure 11).The results showed that, with a similar forest coverage, the CRIRs of the basin-pairs were calculated as 2.48% and 5.08% in the middle subtropical zone, 4.03% and 5.13% in the northern south subtropical zone, 3.35% and 4.15% in southern south subtropical zone, and 2.34% and 3.10% in the north tropical zone, respectively (Figure 12).This shows that there were no significant differences in the influence of climatic zone on canopy rainfall interception in the basins.

Effects of Climate and Forest Cover on Variation Patterns of the CRIR
The forest canopy plays an important role in rainfall interception in the forest ecosystem.However, differentiating the effects of climate and forest cover on the spatiotemporal variation patterns of the CRIR is not well understood for this region.To address this issue, we first selected two basins with similar forest coverage (about 70%) in each of the four climatic zones (e.g., middle subtropical zone, southern south subtropical zone, northern south subtropical zone, and north tropical zone, Table 3) in Guangdong to analyze their variation patterns of CRIR in the same climatic zone (Figure 11).The results showed that, with a similar forest coverage, the CRIRs of the basin-pairs were calculated as 2.48% and 5.08% in the middle subtropical zone, 4.03% and 5.13% in the northern south subtropical zone, 3.35% and 4.15% in southern south subtropical zone, and 2.34% and 3.10% in the north tropical zone, respectively (Figure 12).This shows that there were no significant differences in the influence of climatic zone on canopy rainfall interception in the basins.To eliminate the influence of the climatic zone and evaluate the impact of forest cover change on canopy rainfall interception, we further extracted the forest ecosystem for all integrated basins in each climatic zone, and then analyzed the relationship between forest coverage and annual mean CRIR for each basin.The results showed that the forest coverage in all climatic zones, except for the middle subtropical zone, had a significant linear relationship (p < 0.001) with the CRIR (Figure 13a).The annual mean CRIR increased gradually with an increase in forest coverage in a linear pattern for the southern and northern south subtropical zones with the determination coefficients (R 2 ) of 0.81 and 0.65, respectively.This showed the poor correlation of the middle subtropical and north tropical zones with the R 2 of 0.23 and 0.19, respectively.This means that the forest cover does have a certain impact on canopy rainfall interception, but this effect is not obvious in some areas.Meanwhile, we also analyzed the proportions of different forest types (e.g., broad-leaved forest, coniferous forest, and coniferous and broad-leaved mixed forest) in the same climatic zone (Figure 13b) and found that the broad-leaved forest occupied a higher proportion of the southern and northern south subtropical zones, with a high correlation between the CRIR and forest cover.In contrast, the middle subtropical and north tropical zone had a higher proportion of coniferous forest and coniferous and broad-leaved mixed forest.This indicates that the proportion of broad-leaved forest may play a key role in canopy rainfall interception at the basin scale in Guangdong, but more research is needed in the future.To eliminate the influence of the climatic zone and evaluate the impact of forest cover change on canopy rainfall interception, we further extracted the forest ecosystem for all integrated basins in each climatic zone, and then analyzed the relationship between forest coverage and annual mean CRIR for each basin.The results showed that the forest coverage in all climatic zones, except for the middle subtropical zone, had a significant linear relationship (p < 0.001) with the CRIR (Figure 13a).The annual mean CRIR increased gradually with an increase in forest coverage in a linear pattern for the southern and northern south subtropical zones with the determination coefficients (R 2 ) of 0.81 and 0.65, respectively.This showed the poor correlation of the middle subtropical and north tropical zones with the R 2 of 0.23 and 0.19, respectively.This means that the forest cover does have a certain impact on canopy rainfall interception, but this effect is not obvious in some areas.Meanwhile, we also analyzed the proportions of different forest types (e.g., broad-leaved forest, coniferous forest, and coniferous and broad-leaved mixed forest) in the same climatic zone (Figure 13b) and found that the broad-leaved forest occupied a higher proportion of the southern and northern south subtropical zones, with a high correlation between the CRIR and forest cover.In contrast, the middle subtropical and north tropical zone had a higher proportion of coniferous forest and coniferous and broad-leaved mixed forest.This indicates that the proportion of broad-leaved forest may play a key role in canopy rainfall interception at the basin scale in Guangdong, but more research is needed in the future.

Discussion
Water source conservation is one of the most important vegetation functions in the terrestrial ecosystem.As an interface for atmospheric rainfall in the vegetation ecosystem, vegetation canopy influences the redistribution of rainwater as well as the damping and retention of the water that reaches the soil, affecting runoff dynamics through interception [54,55] .However, previous studies on canopy rainfall interception [55][56][57] have mainly focused on individual plants or groups of plants through in situ observation and are limited by spatiotemporal estimation in larger areas.In this study, we applied a canopy rainfall interception model coupled with MODIS_LAI remote sensing data to understand the spatial and temporal variation patterns of canopy rainfall interception at the basin scale in Guangdong Province, China.
From this study, we found that the CRIR at the basin scale in Guangdong ranged from 0.33% and 11.67%, which is similar to results in other studies in China that ranged from 0% and 12% [58].However, this range is lower than that in forest ecosystems [57,59,60].This is mainly because our study estimated the mean value of the CRIR for all vegetation ecosystems in the basin, which resulted in a lower estimated value when compared with that of a forest ecosystem.It is worth noting that the

Discussion
Water source conservation is one of the most important vegetation functions in the terrestrial ecosystem.As an interface for atmospheric rainfall in the vegetation ecosystem, vegetation canopy influences the redistribution of rainwater as well as the damping and retention of the water that reaches the soil, affecting runoff dynamics through interception [54,55].However, previous studies on canopy rainfall interception [55][56][57] have mainly focused on individual plants or groups of plants through in situ observation and are limited by spatiotemporal estimation in larger areas.In this study, we applied a canopy rainfall interception model coupled with MODIS_LAI remote sensing data to understand the spatial and temporal variation patterns of canopy rainfall interception at the basin scale in Guangdong Province, China.
From this study, we found that the CRIR at the basin scale in Guangdong ranged from 0.33% and 11.67%, which is similar to results in other studies in China that ranged from 0% and 12% [58].
However, this range is lower than that in forest ecosystems [57,59,60].This is mainly because our study estimated the mean value of the CRIR for all vegetation ecosystems in the basin, which resulted in a lower estimated value when compared with that of a forest ecosystem.It is worth noting that the rainfall interception model applied in our study only considered the parameters of the LAI and rainfall, while rainfall interception was also closely related to many other factors such as atmospheric drivers of evaporation (e.g., wind speed, vapor pressure deficit, etc.), canopy characteristics, rainfall intensity, and rainfall duration [61][62][63][64][65].In future research, introducing more meteorological factors into the interception model is necessary to estimate canopy rainfall interception more accurately.
It is important to note that, according to the analysis using the equations of cumulative interception (Equations ( 1) and ( 2)) and CRIR (Equation ( 3)), the spatiotemporal variation of the LAI and rainfall will have a similar influence on the annual and regional variation of vegetation canopy rainfall interception at the basin scale of Guangdong.We also found that the CRIR showed a significant spatial difference pattern that was high in the mountain areas of the eastern and western flanks, but low in the central mountainous region and coastal areas.This spatial pattern was closely related to the LAI of vegetation, which mainly presents a spatial pattern extending from the Pearl River Delta to the periphery, and has little correlation with the spatial pattern characteristics of rainfall that gradually decrease from south to north.Guangdong has many heavy rainfall events, and the leaves tend to reach their maximum saturation of rainfall interception, resulting in vegetation canopy rainfall interception depending mainly on the LAI, but not being restricted by the amount of rainfall.In addition, atmospheric rainfall is the main factor influencing the temporal variation of rainfall interception in the study area.Results showed that the annual mean CRIR in Guangdong continued to decrease from 2004 to 2012, but increased slightly in 2016.However, as time progressed, the difference in CRIRs between the basins decreased.This is because, as the annual mean LAI and annual rainfall increase, the spatial heterogeneity of LAI also increases (Figure 4), while the spatial heterogeneity of rainfall decreases (Figure 2).Combining the effects of the saturation of canopy rainfall interception and rainfall, the spatial homogeneity of annual rainfall magnifies the effects on the spatial pattern of CRIRs, while the spatial heterogeneity of annual mean LAI further diminishes the effects on spatial patterns of CRIR.This has led to decreases in both the CRIR and its spatial variability over time.
Importantly, broad-leaved forests are the main contributor to canopy rainfall interception in Guangdong.By analyzing the CRIRs of the basin-pairs in the four climatic zones with a similar forest coverage, we found that there was no significant difference in the canopy rainfall interception among the four climatic zones, while there was a significant correlation between forest coverage and the canopy rainfall interception within the basin.This means that canopy rainfall interception is far more affected by vegetation coverage than the background climate.In addition, the proportion of broad-leaved forest to the whole area of forest in the basin plays a critical role in the CRIR in Guangdong, which is similar to other results from previous research [55,66,67].However, many studies conducted in various forests have reached the opposite conclusion, demonstrating that the canopy interception capacity in coniferous forests was higher than that in broad-leaved forests [68,69].Researchers have found that the leaf surface area contributes to differences in the interception capacities between broad-leaved forests, coniferous forests, and coniferous and broad-leaved mixed forests [70], and the total annual interception was greater for the species with a denser canopy (e.g., broad-leaved forest).However, since the canopy interception capacity is also determined by other meteorological and canopy parameters, further research is needed on these other factors.
Note also that Guangdong is a typical rainy region in southern China with increased precipitation variability and abnormally frequent rainfall over the past several years, which can be explained at least in part by climate change.Water resource and soil erosion are two of the most crucial ecological concerns in this region and both are closely related to rainfall interception.Therefore, the assessment of vegetation canopy interception is extremely important for characterizing and understanding water cycling, and has scientific and practical implications for land surface cover planning and water and soil conservation.However, the level of economic development and ecological characteristics of various parts of Guangdong Province vary greatly, which requires relevant departments to formulate feasible policy plans according to local conditions.For example, in the part of the Pearl River Delta region that has a high demand for land for construction and economic development, green space system construction should be encouraged.For the mountainous areas in northern Guangdong with relatively large forest areas, it is necessary to strengthen forest conservation in successional stages to improve the quality of forests and form more vegetation coverage.For the Leizhou Peninsula and other coastal areas with serious soil erosion, vegetation types including mangrove forest, broad-leaved forest, and high-density grassland should be the main components of the ecological system in these regions.

Conclusions
In summary, the results of this study provide insights into the spatiotemporal variation in rainfall interception at the basin scale of Guangdong Province, and suggest that forest cover should be increased by adjusting the species composition and increasing the proportion of native broad-leaved species according to local condition in the basin.Furthermore, the results can provide helpful information for water resource management for the local government, and could potentially help assess the impacts of forest ecosystems on water cycling at a larger scale in the context of climate change.
31 N, 109 • 40 -117 • 20 E) and occupies an area of 179,752 km 2 .The climate of Guangdong is tropical and subtropical monsoon and the 12 months can be divided into humid (April to September) and dry (October to March) seasons.The annual mean temperature is 22 • C, and the annual precipitation amount is about 1770 mm, of which 1384 mm (78%) and 386 mm (22%) fall in the wet and dry seasons, respectively [45].

Figure 2 .
Figure 2. Spatial patterns of annual rainfall in Guangdong.Figure 2. Spatial patterns of annual rainfall in Guangdong.

Figure 2 .
Figure 2. Spatial patterns of annual rainfall in Guangdong.Figure 2. Spatial patterns of annual rainfall in Guangdong.

Figure 3 .
Figure 3. Spatial patterns of annual rainfall in the seven major basins in Guangdong.

Figure 3 .
Figure 3. Spatial patterns of annual rainfall in the seven major basins in Guangdong.

Figure 4 .
Figure 4. Spatial patterns of the annual mean leaf area index (LAI) in Guangdong.

Figure 5 .
Figure 5. Spatial patterns of the annual mean LAI in the seven major basins in Guangdong.

Figure 4 .
Figure 4. Spatial patterns of the annual mean leaf area index (LAI) in Guangdong.

Figure 4 .
Figure 4. Spatial patterns of the annual mean leaf area index (LAI) in Guangdong.

Figure 5 .
Figure 5. Spatial patterns of the annual mean LAI in the seven major basins in Guangdong.

Figure 5 .
Figure 5. Spatial patterns of the annual mean LAI in the seven major basins in Guangdong.

Figure 6 .
Figure 6.Spatial pattern of subregion with different rainfall interception in Guangdong in 2004.

Figure 6 .
Figure 6.Spatial pattern of subregion with different rainfall interception in Guangdong in 2004.

Figure 7 .
Figure 7. Variation of monthly mean canopy rainfall interception rates in Guangdong.

Figure 7 .
Figure 7. Variation of monthly mean canopy rainfall interception rates in Guangdong.

Figure 8 .
Figure 8. Spatial patterns of annual mean CRIRs in Guangdong.

Figure 9 .
Figure 9. Spatial patterns of the propensity score (SLOPE) for CRIRs in Guangdong.

Figure 8 .
Figure 8. Spatial patterns of annual mean CRIRs in Guangdong.

Figure 8 .
Figure 8. Spatial patterns of annual mean CRIRs in Guangdong.

Figure 9 .
Figure 9. Spatial patterns of the propensity score (SLOPE) for CRIRs in Guangdong.

Figure 9 .
Figure 9. Spatial patterns of the propensity score (SLOPE) for CRIRs in Guangdong.

Figure 10 .
Figure 10.Spatial patterns of the annual mean CRIR in seven major river basins in Guangdong.

Figure 10 .
Figure 10.Spatial patterns of the annual mean CRIR in seven major river basins in Guangdong.

Figure 11 .
Figure 11.Basin-pairs of the four climatic zones in Guangdong.

Figure 12 .
Figure 12.Annual mean variation of the CRIR in basin-pairs in the four climatic zones in Guangdong.

Figure 11 .
Figure 11.Basin-pairs of the four climatic zones in Guangdong.

Figure 11 .
Figure 11.Basin-pairs of the four climatic zones in Guangdong.

Figure 12 .
Figure 12.Annual mean variation of the CRIR in basin-pairs in the four climatic zones in Guangdong.

Figure 12 .
Figure 12.Annual mean variation of the CRIR in basin-pairs in the four climatic zones in Guangdong.

Figure 13 .
Figure 13.The relationship between the CRIR and forest coverage of all integrated basins in the four climatic zones (a), and the percentage of forest types in different climatic zones (b).

Figure 13 .
Figure 13.The relationship between the CRIR and forest coverage of all integrated basins in the four climatic zones (a), and the percentage of forest types in different climatic zones (b).MSZ: Middle subtropical zone, NSSZ: Northern south subtropical zone, SSSZ: Southern south subtropical zone, NTZ: North tropical zone.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/11/21/2468/sl, Figure S1: Spatial patterns of annual rainfall from 2004 to 2016 in Guangdong, Figure S2: Spatial patterns of annual rainfall from 2004 to 2016 in seven major basins in Guangdong, Figure S3: Spatial patterns of annual mean LAI from 2004 to 2016 in Guangdong, Figure S4: Spatial patterns of annual mean LAI from 2004 to 2016 in seven major basins in Guangdong, Figure S5: Spatial patterns of annual mean canopy rainfall interception rate (CRIR) from 2004 to 2016 in Guangdong, Figure S6: Spatial patterns of annual mean CRIR from 2004 to 2016 in seven major river basins in Guangdong.

Table 1 .
Categorization of canopy rainfall interception.
S v : the canopy rainfall interception (mm).

Table 2 .
Annual canopy rainfall interception (mm) in the different subregions.

Table 2 .
Annual canopy rainfall interception (mm) in the different subregions.

Table 3 .
Areas and percentages of forests in basin-pairs in the four climatic zones.

Table 3 .
Areas and percentages of forests in basin-pairs in the four climatic zones.