Trends in the regional mean time series of Tair, VPD, Tmax, Tmin, Prec, and RH had different characteristics during 1998–2012 (Figure 3). For Tair, the annual mean warming rate during 1998–2012 decreased by −0.14 ◦C/decade, and the decreases in the annual mean Tmax and Tmin, −0.13 ◦C/decade and −0.08 ◦C/decade, respectively, were much weaker and insignificant (Figure 3a,c,d). The annual total precipitation trends show an insignificant decline of −7.66 mm/decade between 1998–2012, which is contrasts markedly with the earlier period (1982–1997) (Figure 3e). The decrease in annual RH, −3.3%/decade, was significant during the warming hiatus (Figure 3f). For the VPD, the significant positive trends of 0.08 hPa/decade were impressive during 1998–2012 (Figure 3b), which is inconsistent with the reduction in other variables during the same time.

**Figure 3.** Regional mean time series of Tair, VPD, Tmax, Tmin, Prec, and RH in Xinjiang during the period 1998–2012 (1982–1997 and 2013–2015 series are shown by gray lines for comparison).

The seasonal SAT trends had different characteristics before and after 1998. The greatest increase in Tair occurred during the spring, 0.52 ◦C/decade, and the summer and autumn had slight warming trends, 0.22 ◦C/decade and 0.14 ◦C/decade, respectively; however, there was a decreasing trend, −1.36 ◦C/decade in winter. This contrasts markedly with 1982–1998, during which there was little change (0.17 ◦C/decade) in summer, while the other three seasons experienced significant warming trends, 0.76 ◦C/decade in spring, 0.93 ◦C/decade in autumn, and 0.98 ◦C/decade in winter (Table 1).

**Table 1.** Change trend for GIMMS3g-NDVI and temperature for the periods 1982–1998 and warming hiatus (1998–2012).

VPD is a direct indicator of atmospheric desiccation strength and is dependent on temperature and RH. The Xinjiang climate experienced remarkable changes during the warming hiatus. The annual mean temperature after 1998 increased by 9.2% compared with that during 1982–1998 (Figure 4). Seasonally, spring had the largest change at 12.4%, followed by autumn at 11.3%, and summer had the weakest change of 3.3% during 1998–2015 relative to 1982–1998 in Xinjiang (Figure 4). For VPD, the relative change was approximately 8.0% before and after 1998, and the seasonal VPD changes showed the greatest increase in spring (12.3%), followed by 6.9% in autumn, while summer showed the smallest change of 6.7% (Figure 4). Thus, on seasonal timescales, spring and autumn became warmer and drier during the warming hiatus period, and summer changed little.

**Figure 4.** Change (%) in temperature during 1998–2015 relative to that in 1982–1998 in Xinjiang. The light gray background shows the spring and autumn.

#### *3.2. Vegetation Browning Trends Based on Remote Sensing Data*

The monthly vegetation NDVI for each grid cell in Xinjiang during 1982–2015 was estimated using the satellite-derived GIMMS3g dataset. To determine the effects of the warming hiatus (1998–2012) on vegetation change, we investigated NDVI trends in Xinjiang before and after 1998. The average vegetation NDVI and temperature trends in Xinjiang before 1998, after 1998, and during the warming hiatus periods, are shown in Figure 5 and Table 1.

**Figure 5.** GIMMS3g-NDVI and temperature trends in Xinjiang from 1982 to 2015. A linear trend was determined for the periods before and after 1998, separately. The light gray background shows the warming hiatus (1998–2012).

The annual NDVI showed a significant greening trend of 0.004/decade (*p* < 0.01) before 1998 but had an insignificant browning trend during 1998–2015 (slope = −0.001/decade, *p* > 0.1), and an NDVI browning rate of −0.002/decade occurred during the warming hiatus period (Figure 5a). Seasonally, the greatest greening rate of 0.007/decade (*p* < 0.01) occurred in summer, followed by a considerable greening rate of 0.06/decade (*p* < 0.01) in autumn, and 0.04/decade (*p* < 0.01) in spring during 1982–1998 (Figure 5b–d, Table 1).

During the warming hiatus period, the reductions in the spring and autumn NDVI, −0.002/decade and −0.003/decade (*p* > 0.05), respectively, were much weaker and insignificant (Figure 5b,d, Table 1). However, the summer NDVI increased significantly, with a greening rate of 0.001/decade (*p* < 0.05) (Figure 4c, Table 1). We also examined the annual and seasonal warming rates (Figure 5). During 1998–2012, the average temperature warming rates reversed at the annual scale (from 0.80 ◦C/decade to −0.14 ◦C/decade) and decreased in spring and autumn (from 0.76 ◦C/decade to 0.52 ◦C/decade and from

0.93 ◦C/decade to 0.14 ◦C/decade, respectively) but continuously increased in summer (from 0.17 ◦C/decade to 0.22 ◦C/decade) (Figure 5, Table 1).

The spatial patterns of the vegetation NDVI were also compared before and during the warming hiatus (Figures 6 and 7). During 1982–1998, approximately 80% of the vegetationcovered area in Xinjiang exhibited a greening trend in spring (Figure 6c,d). In comparison, the greening trend of the spring NDVI reversed in southern Xinjiang, the Junggar Basin, and parts of the Tianshan Mountains, where a significant spring NDVI browning trend occurred during the warming hiatus (Figure 6e,f). In addition, the spring NDVI significantly increased between 1998 and 2012 in the mountainous areas, mainly in the Altay Mountains, northern Tianshan Mountains, and part of the Ili River valley, where there is abundant snowfall in winter.

**Figure 6.** Vegetation NDVI trends (**left**) and M–K significance test (**right**) in spring over Xinjiang, China during 1982–2015 (**a**,**b**) and the two periods of 1982–1998 (**c**,**d**) and 1998–2012 (**e**,**f**).

**Figure 7.** Vegetation NDVI trends (**left**) and M–K significance test (**right**) in autumn over Xinjiang, China during 1982–2015 (**a**,**b**) and the two periods of 1982–1998 (**c**,**d**) and 1998–2012 (**e**,**f**).

In comparison, the NDVI greening trends in autumn were larger than those in spring during 1982–1998, and almost all vegetation-covered areas in Xinjiang exhibited a greening trend in autumn (Figure 7c,d). However, the autumn NDVI showed a significant browning trend during the warming hiatus, and the NDVI browning trend in autumn was greater than that in spring. The browning NDVI area was mainly found in the Tianshan Mountains, northern Xinjiang, and the Kunlun Mountains, which are dominated by natural vegetation (Figure 7e,f). The most dominant vegetation cover type is grassland, which is mainly driven by climate factors. In addition, the autumn NDVI greening trend during 1998–2012 was generally found in the artificial vegetation area, including the northern slope of the Tianshan Mountains and the lower Tarim River Basin.

The lower Tarim River is an ecologically fragile region due to its sensitivities to climate variations and anthropogenic activities. Due to the intensified human influence, a serious decline in vegetation cover occurred in the lower Tarim River during the recent 30 years. The autumn NDVI also showed a significant browning trend in the lower Tarim River during 1982–1998 (Figure 7c,d). The Chinese government has implemented the ecological

water diversion project (EWDP) since 2000 to restore the riparian vegetation in the lower Tarim River. To evaluate the ecological effects of EWDP, Figure 8 shows the time series of the vegetation area and total water conveyance volume over the lower Tarim River from 2000 to 2015 corresponding to the 15 ecological water deliveries. Until the end of 2015, the total water conveyance volume reached 51.11 × 108 m3 (Figure 8). The vegetation processes significantly responded to the water conveyance, the NDVI of the lower Tarim River increased by 28.6%, and the vegetation area increased by 161.2% from 2000 to 2015 (Figure 8).

**Figure 8.** The time series of the vegetation area (km2) and total water conveyance volume (108 m3) over the lower Tarim River during 2000–2015.

#### *3.3. Climatic Controls of Vegetation NDVI*

To examine the effects of climatic factors on vegetation change, we calculated the correlation coefficients between seasonal climatic factors (Tair, Tmax, Tmin, Prec, RH, WS, and SD) and vegetation NDVI in Xinjiang for the 1982–2015, and warming hiatus periods. For the spring vegetation NDVI, we used spring and winter climatic factors, and spring, summer, and autumn climatic factors were used for the autumn vegetation NDVI. The influence of climatic factors on vegetation change varied greatly among the seasonal scales. Generally, the spring NDVI displayed a higher correlation with climatic factors than did the autumn NDVI (Figures 9 and 10).

For the spring vegetation NDVI during 1982–2015, the spring air temperature showed the highest correlation among the examined climatic factors, with correlation coefficients of 0.62 with Tair, 0.61 with Tmax, and 0.56 with Tmin, followed by the spring VPD (correlation coefficient of 0.54, *p* < 0.01) (Figure 9a). During the warming hiatus period, the spring VPD had the highest positive correlation with spring NDVI, followed by spring Tair and Tmax. Moreover, spring precipitation and RH were significantly negatively correlated with spring NDVI compared to any other climatic factor (Figure 9b). The spring browning NDVI area was mainly found in the alpine regions. Spring precipitation is mainly in the form of snowfall in alpine regions of Xinjiang, and snow cover is an essential component of alpine land cover [44]. Snow cover phenology has a strong control on alpine vegetation dynamics, for example, the length of snowing days and snowmelt dates have a significant impact on the NDVI variation in spring over large areas [45]. In summary, spring climatic factors (e.g., air temperature and VPD) were the main controlling factors of the spring NDVI.

Among the various factors, summer precipitation, and RH had the highest correlation with autumn vegetation NDVI during the 1982–2015 and warming hiatus periods, and antecedent and contemporaneous WS were significantly negatively correlated with autumn NDVI during 1982–2015 (Figure 10a,b). The correlation between temperature and NDVI in autumn was much weaker than that in spring. In addition, summer VPD was significantly negatively correlated with autumn NDVI for the warming hiatus period (Figure 10b). Thus, summer climatic factors (e.g., VPD, WS, RH, and precipitation) were significantly correlated with the autumn NDVI.

To illustrate the joint controls of the multiple climatic factors on NDVI change, we used PCA to determine the relationship between vegetation NDVI and climatic factors during the warming hiatus period. Spring temperature (Tair and Tmax) had the highest partial correlation with spring NDVI, while summer precipitation had the highest partial correlation with autumn NDVI. Climatic factors displayed a larger effect on autumn NDVI than spring NDVI. These results were generally like those based on univariate correlation analysis.

**Figure 10.** Correlation coefficients between autumn NDVI and climatic factors. The 95% confidence limits for the correlation are ±0.35 and ±0.51 for 1982–2015 and 1998–2012, respectively. The blue background indicates a positive correlation, and red indicates a negative correlation.

#### *3.4. Impacts of Changes in Drought on Vegetation NDVI*

To investigate the effects of dry and wet conditions on vegetation NDVI changes, we examined the correlation coefficients between vegetation detrended NDVI and the drought index (sc\_PDSI). The vegetation NDVI semi-arid regions tend to respond to drought at longer time scales [46]. For this purpose, we used the sc\_PDSI, which has been reasonably successful at quantifying long-term drought [47], as it uses surface air temperature and a physical water balance model, and can capture the basic effect of global warming on drought through changes in potential evapotranspiration [47]. The procedure followed to calculate the correlations between the standardized vegetation NDVI series in spring and autumn and the series of the sc\_PDSI series, which are illustrated in Figures 11 and 12. For the spring NDVI, we used the spring and winter of the year before (pre-winter) sc\_PDSI, and the summer and autumn sc\_PDSI were used for autumn NDVI. The effects of dry and wet conditions on vegetation NDVI change varied greatly between the spring and autumn. Generally, the autumn NDVI showed a higher correlation with dry and wet conditions than the spring NDVI during 1982–2015 (Figures 11 and 12).

**Figure 11.** Relationships between spring standardized NDVI and standardized sc\_PDSI in Xinjiang. The left graph is a comparison of the time series of spring NDVI and spring sc\_PDSI (**a**), and winter of the year before (pre-winter) sc\_PDSI (**c**) during 1982–2015. The right graph is the scatter plot between the spring NDVI and spring sc\_PDSI (**b**) and between the spring NDVI and previous winter sc\_PDSI (**d**) during 1998–2012. The gray shaded regions in (**a**,**c**) denote the warming hiatus (1998–2012).

**Figure 12.** Relationships between autumn standardized NDVI and standardized sc\_PDSI in Xinjiang. The left graph (**a**) is the time series of autumn NDVI and summer sc\_PDSI and the right graph (**b**) is the scatter plot between the autumn NDVI and summer sc\_PDSI during 1982–2015. The gray shaded regions in (**a**) denote the warming hiatus (1998–2012).

During 1982–2015, the spring NDVI had a slight negative correlation with spring sc\_PDSI (CC = −0.16, *p* > 0.05) and an insignificant positive correlation with pre-winter sc\_PDSI (CC = 0.14, *p* > 0.05) (Figure 11). However, the spring NDVI was significantly negatively correlated with spring sc\_PDSI (CC = −0.71, *p* < 0.01) and slightly negatively correlated with pre-winter sc\_PDSI (CC = −0.40, *p* > 0.05) during the warming hiatus (1998–2012) (Figure 11). In summary, dry and wet conditions in spring and pre-winter were not the key factors controlling spring vegetation changes. For the autumn NDVI during 1982–2015, the summer sc\_PDSI showed a higher correlation (CC = 0.37, *p* < 0.05), followed by autumn sc\_PDSI (CC = 0.33, *p* > 0.05) (Figure 12). However, the correlation between the autumn NDVI and summer (autumn) sc\_PDSI during the warming hiatus was much weaker than that during 1982–2015. Thus, dry and wet conditions in summer were key factors controlling autumn vegetation changes.

#### **4. Discussion**

Based on satellite-derived vegetation records (GIMMS3g-NDVI), a significant greening trend was observed in our study during 1982–1998, and an insignificant browning trend during 1998–2015. Many studies based on satellite records indicated vegetation greening trends in most of the global vegetation areas before the late 1990s [22–24,48], and a recent study showed a pronounced vegetation browning trend over approximately 59% of the global vegetation surface from 1998 to 2015 [22]. Many other studies have reported vegetation greening towards browning in Xinjiang and High Mountain Asia from the late 1990s to the early 2010s [11,12,24,49,50].

Compared to previous studies, we focused on the vegetation change trend in spring and autumn during the recent warming hiatus (1998–2012), which has not been examined in previous studies. According to the definition of the recent warming hiatus, we divided the study period into two shorter periods: before 1998 and after 1998. A previous study using satellite records found that the advancing (or delaying) rate of spring (or autumn) phenology from 2002 to 2012 was faster than that in 1982–1992 [21,51]. Another study using long-term satellite-derived vegetation records and FLUXNET datasets in NH proposed that spring and autumn phenology showed no trends during the global warming hiatus [17]. Our analysis showed that the summer NDVI significantly increased; however, there was an insignificant reduction in the spring and autumn NDVI during the warming hiatus period. We also found that the autumn NDVI browning trend was greater than that in spring. The 56 site FLUXNET database indicated that more records showed significant phenology trends in autumn than in spring [17], and these studies were consistent with our results based on satellite records. We also examined the vegetation change trend in seasons that were consistent with those in seasonal temperature during 1998–2012, and the temperature in spring and autumn decreased but continuously increased in summer.

Our results using the GIMMS NDVI3g database indicated that the vegetation greening trend stalled or reversed during the warming hiatus period. This vegetation browning area is mainly dominated by natural vegetation. In addition, the GIMMS NDVI3g dataset contains uncertainties in data quality, which are caused by noise, sensor degradation, and calibration between different sensors [52].

Many studies have indicated that the spring temperature is the key controlling factor of vegetation change, and rising spring temperatures lead to advancing trends in spring phenology [17,53,54]. Previous studies have also indicated that winter temperature is negatively correlated with spring vegetation NDVI and phenology change [20,55].

Our results suggested that spring NDVI was more largely correlated with spring VPD than with spring temperature during the warming hiatus. Previous studies have shown that an increase in VPD can reduce global vegetation growth [22], and other studies have reported that an increase in VPD is a driver of drought-related forest mortality [56,57]. We also used PCA to confirm that the spring NDVI was mainly dominated by spring temperature, including mean and extreme variables. The VPD represents the difference between saturation vapor pressure (*es*) and the actual vapor pressure (*ea*) at a given temperature [43]. Rising air temperatures may increase the *es*, and the *ea* is constrained at the upper end by *es* so that as air temperature increases, so does the maximum amount of water vapor [58]. Thus, higher VPD is driven by rising air temperature following a nonlinear relationship [58]. Li et al. [43] also indicated that the VPD increase in Xinjiang was jointly determined by the increased saturation vapor pressure from rising temperatures and decreased actual vapor pressure. Thus, rising air temperatures resulted in increased VPD, which corresponded to

the vegetation browning trend. In addition, our results indicated that spring precipitation was significantly negatively correlated with the spring NDVI. Li et al. [29] found that the reversed trends in NDVI after the late 1990s can be attributed to increased warming and reduced precipitation. Yao et al. [12] found that increased drought severity exacerbated the loss of soil moisture and led to a decreased NDVI. These results are consistent with those of our work during the warming hiatus period.

Compared to spring vegetation NDVI, climatic factors displayed a stronger control on autumn NDVI. Wang et al. [17] also proposed that autumn vegetation growth is determined by a combination of climatic factors and phenology statuses. Richardson et al. [59] reported the effects of temperature on autumn vegetation change and neglected other climatic factors [59]. In contrast, our results indicate that summer climatic factors (e.g., VPD, WS, RH, and precipitation) were significantly correlated with autumn NDVI during the warming hiatus. However, the autumn temperature was weakly correlated with the autumn NDVI. Wang et al. [17] found that temperature, VPD, shortwave radiation, precipitation, and soil moisture had strong effects on autumn vegetation among the 56 FLUXNET sites in the NH. These results emphasize the complexity of the impacts on autumn vegetation changes during the warming hiatus period.

Figure 13 shows the mechanisms of the impact of climatic variables on spring and autumn vegetation growth. Our results provide evidence for the different impacts of climatic variables on spring and autumn NDVI during the warming hiatus periods. The difference in spring and autumn temperature has played a key role in recent decades, and a larger increase in air temperature in spring than in autumn has resulted in increased VPD differences in spring and autumn [60,61]. The VPD and extreme drought stress can trigger vegetation growth retardation through hydraulic failure and carbon starvation [57]. In spring, an increased spring VPD enhances evaporative water loss, which triggers stomatal closure. Increased VPD and reduced stomatal conductance result in restricted plant carbon uptake and productivity and drive widespread vegetation browning [62]. Our study suggests the importance of summer climatic variables in controlling autumn vegetation growth in Xinjiang. The autumn NDVI had a higher correlation with dry and wet conditions. Furthermore, increased ET directly influences soil moisture and causes extreme drought. Soil moisture is a direct reservoir for vegetation, and its deficit plays a large role in autumn vegetation browning in Xinjiang. Extreme drought may decrease the amount of water that can be taken up by plant roots, increase catheter embolism, stagnate water transportation, and can drive vegetation browning.

**Figure 13.** Schematic diagram of the impacts of climatic variables on spring and autumn vegetation growth.

The AVHRR-based GIMMS3g-NDVI is a fundamental dataset for the study of global vegetation growth and ecosystem processes. Guay et al. [63] compared the long-term NDVI trends in high northern latitudes (>50◦ N), between the GIMMS3g dataset and four commonly used NDVI datasets. and found that the marked differences were both in mean seasonal NDVI and long-term trends. The vegetation greening was more prominent in GIMMS 3g-NDVI than in the older dataset, especially in the tundra, needle-leaved evergreen and mixed-leaf forests [63–65], while the browning was across sizeable areas of boreal forest [66,67]. Consistent with the MODIS record, the spatial patterns of GIMMS 3g-NDVI trends showed roughly half greening and half browning across high northern latitudes [63]. In the High Mountain Asia region, the GIMMS 3g-NDVI showed a browning trend, whereas MODIS and SPOT exhibited a greening trend during 2001–2015 [68]. Thus, it was necessary to provide a comprehensive comparison with different NDVI datasets, as well as to validate these NDVI data against in situ measurements of ecosystem processes.

#### **5. Conclusions**

This study focused on the latest changes in spring and autumn vegetation growth in Xinjiang and their climate-driving mechanisms, using satellite-derived vegetation records (GIMMS3g-NDVI) as indicators. The results indicated that the vegetation NDVI showed a significant greening trend during 1982–1998 and an insignificant browning trend during 1998–2015. We focused on the vegetation change trends in spring and autumn during the recent warming hiatus (1998–2012), which had not been examined in previous studies.

During the recent warming hiatus, the summer NDVI increased significantly; however, there were insignificant reductions in the spring and autumn NDVI. The autumn NDVI browning trend in Xinjiang was greater than that in spring. The seasonal vegetation change trend was consistent with the seasonal temperature, and the temperature slowed down in spring and autumn but continuously increased in summer during the warming hiatus.

The correlation coefficients between climatic factors and vegetation NDVI showed that the climatic factors controlling vegetation change varied greatly at the seasonal scale. Spring NDVI displayed a stronger correlation with climatic factors than did the autumn NDVI. During the warming hiatus, the spring VPD had the highest positive correlation with the spring NDVI, followed by spring temperature. In addition, spring precipitation has a significantly negative correlation with the spring NDVI. Thus, spring climatic factors (e.g., air temperature and VPD) were the main controlling factors of the spring NDVI.

Compared to the spring vegetation NDVI, climatic factors displayed a stronger control over the autumn NDVI. Summer climatic factors (e.g., VPD, WS, RH, and precipitation) were significantly correlated with the autumn NDVI during the warming hiatus. However, the autumn temperature was weakly correlated with the autumn NDVI. These results emphasize the complexity of the impacts on autumn vegetation changes during the warming hiatus.

We also provided evidence for the different impacts of climatic variables on spring and autumn vegetation NDVI during the warming hiatus period. The larger increase in air temperature in spring than in autumn resulted in increased VPD differences in spring and autumn. Increased spring VPD can trigger spring vegetation browning through mechanisms of carbon starvation; however, summer drought may decrease the amount of water that can be taken up by plant roots and can drive autumn vegetation browning.

**Author Contributions:** Conceptualization, M.L., J.Y. and J.Z.; methodology, M.L.; formal analysis, M.L. and J.G.; writing—original draft preparation, M.L.; writing—review and editing, M.L., J.Y., J.G. and J.Z.; funding acquisition, J.Y. and J.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the National Natural Science Foundation of China (U1903113), the Major Program of National Social Science Foundation of China (No. 17ZDA064), and the Xinjiang Tianshan Cedar Project (No. 2020XS04).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author or the reference website.

**Acknowledgments:** We thank "anonymous" reviewers for their original insights. We are also immensely grateful to the editor for their comments on the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Spatial-Temporal Variation in Paddy Evapotranspiration in Subtropical Climate Regions Based on the SEBAL Model: A Case Study of the Ganfu Plain Irrigation System, Southern China**

**Guangfei Wei 1, Jingjing Cao 2, Hua Xie 1, Hengwang Xie 3, Yang Yang 4, Conglin Wu 4, Yuanlai Cui 1 and Yufeng Luo 1,\***

**Abstract:** The surface energy balance algorithm for land (SEBAL) is a commonly used method for estimating evapotranspiration (ET) at a regional scale; however, the cloudy and rainy characteristics of subtropical monsoon regions pose a greater challenge for estimating paddy field ET based on remote sensing technology.  Subsequently, we evaluated the applicability of the SEBAL model for estimating the ET of paddy fields at the daily scale; derived the interannual variation (2000–2017) characteristics of early, middle, and late rice ET; and finally analyzed the spatial distribution patterns of rice in different hydrological years. The results demonstrated that: (1) the SEBAL model estimated ET accurately on a daily scale, with *R*2, NSE, and RMSE values of 0.85, 0.81, and 0.84 mm/day, respectively; (2) the ET of paddy fields in the irrigated area was higher in July and August and the interannual trend of ET of early rice was not obvious, with a declining trend observed in middle rice and late rice from 2000 to 2009, which was followed by an increasing trend from 2009 to 2017; and (3) variations in the spatial distribution of ET were significant for early and late rice at different precipitation levels and less obvious for middle rice in wet years but significant in dry years. Overall, this study verified the applicability of the SEBAL model for estimating ET in paddy fields in subtropical regions and provided a basis and reference for the rational allocation of water resources at a regional scale.

**Keywords:** evapotranspiration; paddy rice; SEBAL; interannual variation; spatial distribution; Ganfu Plain irrigation system

#### **1. Introduction**

Paddy rice is a dominant staple food that feeds over 50% of the population worldwide. However, paddy rice needs a large amount of water to maintain high production. Statistics have shown that paddy field irrigation accounts for 80% of agricultural water consumption in southern China [1]. Evapotranspiration (ET) represents the sum of vegetation transpiration and soil evaporation, which is the major component of farmland water consumption, and it includes the amount of evaporation on the ground surface and plant surface and the amount of transpiration that occurs at the plant surface and inside the plant. ET not only plays a huge role in hydrological cycles and energy balance systems, but also represents a crucial link between ecological and hydrological processes; furthermore, it is a significant

**Citation:** Wei, G.; Cao, J.; Xie, H.; Xie, H.; Yang, Y.; Wu, C.; Cui, Y.; Luo, Y. Spatial-Temporal Variation in Paddy Evapotranspiration in Subtropical Climate Regions Based on the SEBAL Model: A Case Study of the Ganfu Plain Irrigation System, Southern China. *Remote Sens.* **2022**, *14*, 1201. https://doi.org/10.3390/rs14051201

Academic Editors: Jingzhe Wang, Zhongwen Hu, Yangyi Wu and Jie Zhang

Received: 26 January 2022 Accepted: 25 February 2022 Published: 28 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

parameter for calculating the regional water balance and a convenient means to retrieve water amounts applied for irrigation over anthropized regions [2–4]. Water consumption dominated by ET is the final form of water resource consumption, which often approximates the regional water resource consumption [5]. With economic and social development and population growth, the problem of water resource shortages is becoming increasingly prominent. The unbalanced spatial and temporal distribution of water resources is one of the challenges of water resource management. Therefore, understanding the spatial and temporal variability in ET is of great significance for rationally allocating regional water resources and scheduling irrigation [6].

With the continuous development of remote sensing (RS) technology in recent decades [7–9], satellite sensors with multiple spatial, temporal, and spectral resolutions (e.g., Land Remote Sensing Satellite (Landsat), Moderate Resolution Imaging Spectroradiometer (MODIS), Advanced Very-High-Resolution Radiometer (AVHRR), and Sentinel) have become ideal tools for estimating ET at a regional scale [10–13]. Various models have been reported for estimating ET with RS datasets [14–16]. According to different theoretical bases, RS-based ET models can be roughly classified as: (i) empirical statistical methods, which have the advantages of requiring less data and being easy to use; however, they have significant regional limitations because their statistical parameters need to be recalibrated when used in different regions [17,18]; (ii) feature space methods, which refers to the method of obtaining feature points and interpolating the ET of the feature points to derive the ET of each image element based on each surface parameter obtained from remote sensing images—such as surface albedo, vegetation index, surface temperature, etc.—presenting a certain shape on the scatter plot, and using that shape to obtain the ET of the feature points [18,19]; (iii) residual surface energy balance methods, which calculate the net radiation, soil heat flux, and sensible heat flux of the surface energy balance equation and finally calculate the latent heat flux as the residual term, which can be divided into single-source models and dual-source models according to whether soil evaporation and vegetation transpiration are distinguished. Studies have shown that vegetation cover can affect the water cycle and accurate estimation of ET [20,21]. Since the dual-source model considers the parameters of vegetation and soil separately, it can describe the surface fluxes of nonuniform subsurfaces in greater detail; however, it is difficult to obtain the corresponding parameters for the model [22]. The uncertainty and errors of the dual-source model will increase, and the estimation results may not be more accurate than those of the single-source model [23]. In contrast, the single-source model simplifies the treatment of the subsurface, which makes the model parameters easier to obtain and more flexible to apply.

As a typical single-source model, the surface energy balance algorithm for land (SE-BAL) based on the surface energy balance has been applied to estimate ET in different regions with multiple data sources [24–30], since it does not excessively rely on groundassisted data. The accuracy of ET estimated by the SEBAL model can reach 85% for the daily ET, 95% for the whole-season ET at a field scale, and 96% for the annual ET in large watersheds [31]. The SEBAL-based ET model has been gradually applied to water use and water management in irrigation districts in recent years. The SEBAL model has been used to estimate ET across various land uses in Nansi Lake Wetland of China using Landsat-7 ETM+ images, and the results showed that the SEBAL model can calculate the spatial distribution of ET despite little ground-based weather data over large areas [32]. Gao et al. (2020) [33] estimated the actual ET of the Loess Plateau in China from 1990 to 2015 using the MODIS and AVHRR datasets based on the SEBAL model, and the *R*2 values were above 0.87 in 161 subbasins, which indicated the simulation results are reliable and effective; meanwhile, they obtained an interesting finding that the SEBAL model overestimates regional ET, the main reason mainly being due to the overestimation of regional latent heat or underestimation of sensible heat. Ding (2020) [34] compared the performance of the SEBAL model and Surface Energy Balance System model (SEBS, a dual-source model) in estimating regional ET, and results showed that the SEBAL model is suitable for estimating

ET in plain areas, while the SEBS model is suitable in composite terrain. The above studies show that the SEBAL model has been used for estimating ET at field scale to regional scale with encouraging results, but the application of the SEBAL model to rice-dominated irrigation areas has not been sufficiently investigated.

The rice-planting areas are mainly distributed in the monsoon regions of East Asia, Southeast Asia, and South Asia, as well as the tropical rainforests of Southeast Asia. The SEBAL model has been applied to several regions in China: the northeast Sanjiang River Basin [35,36], the North China Plain [37], the Loess Plateau [33,38], the agro-pastoral ecotone of Northwestern China [39], the Taklimakan Desert [40], and the Yellow River Delta [21]. Additionally, the SEBAL model has been widely applied over a multitude of different climatic conditions, including subtropical areas [41–43], but the cloudy and rainy weather in the subtropical climate zone—which is represented by the southern region of China, where mixed double-season early and late rice and single-season middle rice are cultivated—presents a greater challenge for accurate estimations of ET using the SEBAL model.

Therefore, the Ganfu Plain irrigation system (GFPIS) located in southern China, which is a subtropical monsoon climate region, was selected as the study area. The objectives of this study are to: (1) assess the applicability of the SEBAL model for estimating ET in rice-dominated irrigation areas; (2) analyze the daily variation in paddy rice ET over the growing season; and (3) analyze the interannual variation (2000–2017) characteristics of early, middle, and late rice ET and the spatial distribution pattern in different hydrological years. Hence, this study estimated ET of the rice-dominated irrigation areas in southern China using Landsat and the SEBAL model, and derived the spatial-temporal variation patterns of ET in the region, which will provide a reference for regional farmland water consumption and a decision-making basis for the rational and efficient allocation of regional water resources.

#### **2. Material and Methods**

#### *2.1. Study Area*

The GFPIS (8◦10 –28◦48 N, 115◦48 –116◦17 E) is the largest diversion irrigated area south of the Yangtze River in China (Figure 1) and one of the main grain-producing areas in Jiangxi Province, with an effective irrigated area of approximately 67,000 Hm2. The main crop planted in this region is rice, which can be divided into double-season early and late rice and single-season middle rice. The growing season is shown in Figure 2 [44]. Among the rice crops, double-season early-late rice may only be planted with early-late rice, which is rotated with cash crops, such as vegetables, melons, and fruits.

The distribution of water conservancy projects in the GFPIS is shown in Figure 1d, with a total of 10 major trunk canals. With the expansion and economic development of Nanchang city in recent years, the water consumption of urban residents and industry in irrigated areas has been growing and the demand for water resources from various water-using units has been increasing; therefore, the contradiction between water supply and consumption in irrigated areas will be more prominent.

The GFPIS is located in a humid subtropical climate zone, with an average annual temperature of 17–18 ◦C, an average annual wind speed of 3.4 m/s, an average annual sunshine of 1603.4 h, a frost-free period of 260–280 days, an average annual precipitation of approximately 1662.5 mm, and an average annual land evaporation of approximately 800 mm. The irrigated area has four distinct seasons, with low temperatures and rain in spring, greater precipitation and flooding in late spring and early summer, high temperatures and little rain in summer, mild temperatures but little rain in autumn, and cold temperatures in winter.

**Figure 1.** Overview of study area. (**a**) Geographical location of Jiangxi Province; (**b**) Location of the study area; (**c**) DEM of the study area; and (**d**) Distribution of the water conservancy project in the GFPIS.

**Figure 2.** Common growth stages of early, late, and middle rice in the GFPIS: 1–sowing; 2–seeding; 3–transplanting; 4–tillering; 5–booting; 6–heading; 7–milk; 8–mature; and 9–harvest.

#### *2.2. Data Used*

The input data (Table 1) in the SEBAL model include satellite imagery, meteorological data, and digital elevation data (DEM).

**Table 1.** Information of datasets used.

The Landsat images have the same spatial resolution and continuous temporal coverage. All of the Landsat images are level 1 T data, with terrain corrections referencing the World Reference System-2 and intercalibrated across the various Landsat sensors; thus, these images are suitable for time series processing analysis. According to the rice growing season in the GFPIS, satellite imagery acquisition occurred between 20 April and 31 October from 2000 to 2017. All cloud-free Landsat 5 TM, 7 ETM+, and 8 OLI images (Path 121 and Row 40) from the U.S. Geological Survey (USGS) Earth Explorer data portal

(http://earthexplorer.usgs.gov/, accessed on 10 January 2022) were downloaded. A total of 63 images were used, and Figure 3 shows the temporal distribution of the available Landsat images.

**Figure 3.** Temporal distribution of the available Landsat images for the rice growing season from 2000 to 2017.

The meteorological data in this study are from the China Meteorological Data Sharing Service Network (http://data.cma.cn, accessed on 10 January 2022), and the Daily Climatic Dataset (V3.0) of the China International Ground Exchange station was selected as the basic meteorological data. Meteorological data were collected from Nanchang Station from 2000 to 2017, and they included the maximum temperature, minimum temperature, average temperature, average wind speed, relative humidity, and sunshine duration.

DEM data are derived from SRTM data, which are jointly measured by NASA and the National Geospatial-Intelligence Agency (https://earthexplorer.usgs.gov/, accessed on 10 January 2022). The spatial resolutions of SRTM data are 30 and 90 m. SRTM data from radar images have the highest resolution and represent the most widely used terrain data thus far; furthermore, global SRTM data have a unified coordinate system. In this study, DEM data with 30 m resolution were mosaicked together with two sheets and then cropped to a rectangle containing the range size using the GFPIS vector file, as shown in Figure 1c.

#### *2.3. Selection of Typical Hydrological Years*

Precipitation and irrigation are the main sources of water in paddy fields, and ET is the main consumer of water in paddy fields; thus, it is also the basis of irrigation management. Therefore, it is of great guiding significance for water resource management to investigate the distribution characteristics of ET under different precipitation conditions in irrigated areas.

First, the total precipitation during the growing seasons of early, middle, and late rice from 2000 to 2017 was statistically analyzed, and the results are shown in Figure 4. Overall, the total precipitation during the growing season of early rice was the highest and ranged from 351.8 mm to 1025.9 mm, with the minimum value occurring in 2009 and the maximum value in 2010. The second highest was the total precipitation during the growing season of middle rice, with a minimum value of 216.2 mm (2016) and a maximum value of 563.8 mm (2014). The total precipitation of late rice was the lowest and ranged from 58.9 mm to 401.2 mm. These results indicate that precipitation is abundant in the early rice growing season and less abundant in the late rice growing season, and the total precipitation in the middle rice growing season is between that of early rice and late rice.

**Figure 4.** Interannual variation in precipitation during the early, middle, and late rice growing seasons from 2000 to 2017.

In this study, the precipitation data for the GFPIS from 1952 to 2017 were sorted by frequency, the Pearson type III curve was used to fit the line, and the precipitation values corresponding to each frequency case were obtained by fitting the line. Generally, the precipitation values corresponding to P = 25%, 50%, and 75% can be taken as the design values for wet, normal, and dry years, respectively, and then the real years equal to or close to the design values can be selected as the corresponding wet, normal, and dry representative years in the existing actual measurement series. The wet years selected in this paper were 2002 and 2014; the normal years were 2004, 2008, and 2010; and the dry years were 2007, 2008, and 2016.

Consistent with the previous conclusion, the total precipitation during the growing season of rice in different hydrological years was ranked as early rice, middle rice, and late rice. The precipitation during the growing season of early rice was higher than 400 mm in different hydrological years; the precipitation during the growing season of middle rice was as high as 564 mm in wet years, 365 mm in normal years, and only 216 mm in dry years; and the precipitation during the growing season of late rice was 401 mm in wet years, 255 mm in normal years, and as low as 136 mm in dry years.

#### *2.4. SEBAL Algorithm Description*

Since the SEBAL model was proposed by Bastiaanssen, it has been successfully applied to a variety of ecosystems through continuous improvement and development [45]. The SEBAL model calculates the latent heat flux by estimating the net radiation, soil heat flux, and sensible heat flux to obtain the instantaneous ET. The theoretical basis of the SEBAL model is the surface energy balance formula, and the calculation formula is

$$
\lambda \mathbf{ET} = \mathbf{R\_n} - \mathbf{G} - \mathbf{H} \tag{1}
$$

where λET is the latent heat flux (W/m2), *Rn* is the net radiation flux (W/m2), *G* is the soil heat flux (W/m2), and *H* is the sensible heat flux (W/m2).

*Rn* is the main energy source in the transport and exchange of energy, momentum, and water on Earth's surface, and it is calculated as follows:

$$R\_{\mathbf{n}} = (1 - \alpha\_0)\mathbf{R}\_{\mathbf{s}} + \mathbf{R}\_{\mathbf{L}} \downarrow - \mathbf{R}\_{\mathbf{L}} \uparrow - (1 - \varepsilon\_0)\mathbf{R}\_{\mathbf{L}} \downarrow \tag{2}$$

$$R\_{L\_{\downarrow}} = \varepsilon\_a \times \sigma \times T\_a{}^4 \tag{3}$$

$$R\_{L\uparrow} = \varepsilon\_0 \times \sigma \times T\_s^4 \tag{4}$$

where *Rs* is the incoming shortwave radiation (W/m2), *RL* ↓ is the incoming longwave radiation, *RL* ↑ is the outgoing longwave radiation (W/m2), (1 − *ε*0)*RL* ↓ is the incident longwave radiation emitted through the surface (W/m2), *α*0 is the surface albedo (dimensionless), *εa* is the atmospheric radiance, *ε*0 is the thermal surface emissivity (dimensionless), *Ta* is the temperature of the air near the surface, and *Ts* is the surface temperature (K).

$$
\varepsilon\_0 = 1.009 + 0.047 \ln(\text{NDVI}), \text{ NDVI} \ge 0.002 \tag{5}
$$

$$
\varepsilon\_0 = 1, \cdot NDVI < 0.002,\tag{6}
$$

$$T\_s = \mathcal{K}\_2 / \ln\left(\frac{\varepsilon\_{NB} \* K\_1}{R\_c} + 1\right) \tag{7}$$

where *K*1 and K2 are the constant values (K), *Rc* is the corrected thermal radiance of the surface, and *εNB* is the narrow band surface emissivity.

$$NDVI = \frac{\rho\_{NIR} - \rho\_{RED}}{\rho\_{NIR} + \rho\_{RED}} \tag{8}$$

where *ρNIR* and *ρRED* are DN values of the near-infrared band and red band, respectively.

It should be pointed out that the radiance changes at different times in the year and between years as a result of different solar incidences over the surface, especially in rough terrain. The ratio properties of NDVI enable the great mass of noise due to changing sun angles, clouds, topography, or shadow to be cancelled, making NDVI less susceptible to illumination conditions [46].

$$\mathfrak{a}\_0 = \sum \gamma\_i \mathfrak{a}\_i \tag{9}$$

where *αi* and *γi* are the surface reflectance and the weight coefficient of band *i*, respectively. *G* is the soil heat flux, and its equation is

$$G/R\_n = T\_s \ast \left( 0.0038 + 0.0074 \ast a \right) \ast \left( 1 - 0.98 \ast NDVI^4 \right) \tag{10}$$

where *Ts* represents the surface temperature (C); *α* represents albedo; and *NDV I* is the normalized difference vegetation index (dimensionless).

$$H = \frac{\rho\_{\text{air}} \ast \mathbb{C}\_p \ast dT}{\gamma\_{\text{air}}} \tag{11}$$

where *ρ*air is the air density (kg/m3), *Cp* is the air specific heat (J/kg\*K), *γah* is the aerodynamic resistance against heat transfer (s/m), and *dT* is the temperature difference (T1 − T2) at heights of Z1 and Z2 ( ◦C).

$$\gamma\_{\text{ah}} = \frac{\ln(\frac{z\_2}{z\_1})}{k \ast \mathcal{U}\_\ast} \tag{12}$$

where *k* is the von Karman constant (0.41) and *U*∗ is the fraction velocity at heights of 2 m and 0.1 m.

$$\mathcal{U}\_{\*} = \frac{k \ast \mathcal{U}\_{200}}{\ln(\frac{z\_{200}}{z\_{0\text{uv}}})} \tag{13}$$

where *z*0*m* is the surface momentum roughness height (m). The equation for *z*0*m* is

$$Z\_{0m} = \exp(\frac{a \ast N DVI}{a} + b) \tag{14}$$

where *a* and *b* are empirical coefficients.

$$dT = aT\_s + b \tag{15}$$

where *a* and *b* are linear parameters, and *dT* is the temperature difference (T1 − T2) at heights of Z1 and Z2 ( ◦C), with Z1 = 0.1 m, Z2 = 2 m usually used in the SEBAL model. The purpose of defining the zero-plane displacement is to reduce the obstruction of wind by ground vegetation.

In Equation (1), ET represents instantaneous evapotranspiration. The evaporative fraction (Λ) was taken to extrapolate the instantaneous ET to the daily ET in the SEBAL model, and the equations are

$$
\Lambda = \frac{ET}{R\_{\text{fl}} - \text{G}} \tag{16}
$$

$$ET\_{24} = \frac{\Lambda \cdot (R\_{24} - G\_{24})}{\lambda} \tag{17}$$

where *λ* is the latent heat of vaporization (MJ kg−1); ET, *Rn*, and *G* are the instantaneous latent heat flux, net radiation and ground heat flux in each pixel, respectively; *R*24 is the daily net radiation (W m−2); and *G*24 is the daily soil heat flux (W m−2).

#### *2.5. Calculation of ET in Paddy Fields*

The paddy rice fields of the GFPIS were extracted in our previous study [44], and we took the average ET value calculated based on the extracted rice image elements as the daily average paddy rice ET value calculated in Section 3.2.

For the days without accessible RS images, the daily ET can be estimated by linearly interpolating over periods between two consecutive images [47,48]. The specific steps are as follows: (1) Calculate the daily crop reference ET (ET0) using the Penman–Monteith method recommended by the FAO [49]; (2) Calculate the ratio of ET24 to ET0 in two available RS images; and (3) Linearly interpolate the two-spoke image to determine the ratio on the missing dates. The ET24 of the missing date can be obtained by multiplying the ratio by the ET0 on the missing date [50]. The Penman model was modified by the FAO in 1998, and the modified Penman model has a solid theoretical basis and high universality; thus, it is currently the most accurate method of calculating ET0. The calculation formula of the Penman model is as follows:

$$ET\_{0, \, PM} = \frac{0.408\Delta (R\_n - G) + \gamma \frac{900}{T + 273} \mu\_2 (e\_s - e\_a)}{\Delta + \gamma (1 + 0.34 \mu\_2)}\tag{18}$$

$$ET\_{\text{scason}} = \sum\_{i=\text{start}}^{i=\text{end}} ET\_{\text{daily\\_i}} = \sum\_{i=\text{start}}^{i=\text{end}} ET\_{\text{i}} \times ET\_{\text{0\\_PM\\_i}} \tag{19}$$

#### *2.6. Evaluation of Model Accuracy*

In this study, the daily average ET observed based on the eddy covariance of the Jiangxi Province Irrigation Experiment Station (116◦00 03 E, 28◦26 27 N) in the GFPIS from 2016–2017 was compared with the results calculated by the SEBAL model. To evaluate the prediction accuracy of the model for ET, the coefficient of determination (*R*2), root-meansquare error (RMSE), and Nash–Sutcliffe efficiency coefficient (NSE) were chosen to analyze the ET derived by the model, and the calculation process is shown in Equations (20)–(22). A reliable model usually has high *R*2 and NSE values and low RMSE values [51,52].

$$R^2 = 1 - \frac{\sum\_{i=1}^{n} \left( ET\_{\text{EC\\_i}} - ET\_{\text{sat\\_i}} \right)^2}{\sum\_{i=1}^{n} \left( ET\_{\text{EC\\_i}} \right)^2} \tag{20}$$

$$RMSE = \sqrt{\frac{\sum\_{i=1}^{n} \left( ET\_{\text{EC\\_j}} - ET\_{\text{sat\\_j}} \right)^2}{n}} \tag{21}$$

$$NSE = 1 - \left[\frac{\sum\_{i=1}^{n} \left(ET\_{EC\\_i} - ET\_{sat\\_i}\right)^2}{\sum\_{i=1}^{n} \left(ET\_{EC\\_i} - ET\_{EC\\_i}\right)^2}\right] \tag{22}$$

where *ETEC*\_*i* is the *ET* observed by vorticity on day *i*, *ETsat*\_*i* is the *ET* estimated by remote sensing on day *i*, *ETEC*\_*i* is the mean *ET* value observed by vorticity, and *n* is the total number of observed samples.

In this study, ArcGIS (ESRI, West Redlands, CA, USA) and ENVI (Environment for Visualizing Images, Exelis Visual Information Solutions, Boulder, CO, USA) were used to process the RS images and PyCharm (JetBrains, Prague, Czech Republic) was used to construct the SEBAL model using the Python programming language.

#### **3. Results and Discussion**

#### *3.1. Estimation Accuracy of the SEBAL Model for Daily ET in the GFPIS*

In this paper, Landsat images on typical days with cloudless skies or less cloud cover were selected to analyze the estimated surface ET results in the growing season of 2017, which were 25 May, 12 July, 14 September, and 30 September 2017. Due to cloud cover and other influences, the processing results of some pixels were still null, as shown in Figure 5.

**Figure 5.** Daily ET estimated in the GFPIS in 2017 by the SEBAL model.

In terms of time series, the maximum ET value in the GFPIS was 7.6 mm/day on July 12, followed by those on 25 May and 14 September, and the minimum daily ET was 5.6 mm/day on 30 September. The ET of the farmland area on 12 July and 14 September was higher, and the ET of the farmland area on 30 September was lower than that on 14 September. In terms of the four seasons, the maximum and minimum ET values were observed in summer and winter, respectively; moderate values were observed in spring and autumn; and higher values were observed in spring than in autumn.

The ET was estimated by the SEBAL model in 2016 and 2017 and the ETs observed based on the eddy covariance at the experimental station are shown in Figure 6. The RMSE of the comparison between the model estimation results and the observation results was 0.84 mm/d, the *R*2 was 0.85, and the NSE was 0.81. Figure 6 shows that the estimated ET results in 2016 and 2017 were consistent with the observed eddy covariance values. In general, the ET estimated by the SEBAL model was of high accuracy, and the adaptability and stability of the SEBAL model in the irrigated area of GFPIS were satisfactory.

**Figure 6.** Validation of the estimated ET.

The results of this study indicate that the SEBAL model has reasonable accuracy when predicting the ET of subtropical regions. In Tan et al.'s previous study in 2021 [42], the SEBAL model was applied to estimate the ET characteristics in a subtropical region (Huaihe River Basin, China) and the SEBAL model performance was evaluated by fitting the regression with the daily reference ET calculated by multiple theoretical methods. The results showed that the bias between the ET estimated by the SEBAL model and daily reference ET was less than 1.5%. The results of this study have a slightly higher error than that in the previous study. The bias of the SEBAL model's estimation of ET is mainly due to the following factors: (1) The SEBAL model relies on instantaneous surface temperature, and the selection of hot and cold pixels in the imagery is subjective, which has a great impact on the accuracy of ET estimation; (2) The SEBAL model has been proven to be applicable for estimating ET in flat surfaces with relative accuracy, while it tends to have a high bias when applied to mountainous areas [53], which may lead to some increase in bias due to the presence of some hilly areas in our study area; and (3) Studies have shown that even a thin layer of clouds can lead to large errors when calculating sensible heat fluxes [54].

Landsat images with a resolution of 30 m can reflect spatial differences in the ET distribution. Moreover, ET is greatly affected by different ground object types. The greatest evapotranspiration occurred from water bodies—mainly lakes and rivers—whose ET was approximately 5 mm/day. The daily ET of the urban construction area and beach sand was the lowest at only 0–3 mm/day. The ET of forestland and farmland varied with different growth stages. Thus, the estimated spatial distribution of evapotranspiration in the irrigated area is consistent with the actual evapotranspiration.

#### *3.2. Variation in the Daily Rice ET during the Growing Season*

The daily variation characteristics of ET from paddy fields in the GFPIS were analyzed based on the ET of double-season rice in 2016 and 2017, and the results are shown in Table 2.

**Table 2.** SEBAL-estimated daily average evapotranspiration from early and late rice in 2016 and 2017.

In June 2016, as the rice grew and changed from the booting stage to the heading stage. The average ET value increased from 4.69 mm/day on June 7 to 7.43 mm/day on 23 June. On 9 July, as the double-season early rice was in the mature stage, the ET was weak and had a value of 6 mm/day. On 25 July, as the double-season late rice was about to enter the tiller stage, the ET increased to 7.08 mm/day. In September, the ET of double-season late rice was approximately 5.3 mm/day; and in October, it decreased to 4.15 mm/d as double-season late rice entered the mature stage.

In May 2017, as rice growth and transpiration increased, the daily ET value of rice also increased gradually from 5.87 mm/day on 9 May to 6.31 mm/day on 25 May. In July and August, the peak of the double-season early and late rice growing season occurred, the temperature in the GFPIS was high, and the daily ET of rice fields also maintained a high value, which was higher than 6 mm/day. On 14 September, when the late rice was in the heading stage, transpiration remained at 6.09 mm/day. However, the transpiration intensity gradually decreased on 30 September when the rice was in the milking stage and water was lacking in the paddy field due to evaporation by the sun; therefore, the transpiration value decreased significantly to 4.85 mm/day. In November, the transpiration value reached a minimum of 2.69 mm/day when there was no vegetation in the paddy field.

The results showed that the seasonal ET exhibited significant changes and was closely related to the climate. Overall, excluding November and December after rice harvest, the regional daily average double-season early and late rice ET values in 2016 and 2017 showed a pattern higher than 6 mm/day in July and August and lower than 5 mm/day in late September and October, while the daily average paddy ET values varied in the range of 5 to 6 mm/day for the rest of the year.

#### *3.3. Interannual Variation in ET and Precipitation in the Rice Growing Season*

Due to cloudy and rainy weather during the growing season of early rice, limited remote sensing image data are available, and no data are available to calculate the ET of the growing season of early rice in 2003, 2008, 2010, and 2015. Due to the lack of data, we cannot obtain the distribution map of early, middle, and late rice in each year. According to our extraction results for paddy rice in the GFPIS in a previous study [37], the rice planting structure changed, with 2010 as the boundary; according to the available quantity and image quality of satellite data sources during the rice growing period, the middle and double-season rice maps of 2003 were used to calculate the ET of paddy fields from 2000 to 2009 and the middle and double-season rice maps of 2017 were used to calculate the ET of paddy rice from 2010 to 2017. The division of the rice growing season was determined by the observation data from the Jiangxi Irrigation Experiment Station, where the early rice growing season was from 22 April to 10 July, the middle rice growing season was from 26 June to 7 October, and the late rice growing season was from 19 July to 24 October.

Figure 7 shows the interannual variation in ET and precipitation in the rice growing season by the SEBAL model combined with the paddy rice distribution map. The interannual variation in evapotranspiration in the GFPIS from 2000 to 2017 shows that the ET of middle rice is generally the largest, with a multiyear average value of 434 mm; followed by the ET of the late rice growing season, with a multiyear average value of 380 mm; and the ET of early rice is the smallest, with a multiyear average value of 287 mm. The ET in the early rice growing season ranged from 263 mm to 349 mm, with maximum and minimum values occurring in 2013 and 2014, respectively. The average ET value during the growing season of middle rice ranged from 378 mm to 520 mm, with maximum and minimum values occurring in 2016 and 2009, respectively. The mean ET value during the growing season of late rice ranged from 324 mm to 457 mm, with maximum and minimum values occurring in 2016 and 2012, respectively.

**Figure 7.** Interannual variation in ET and precipitation in the rice growing season in the GFPIS.

The interannual variation trend of ET during the rice growing season of early rice was not obvious, but the variation trend of ET of middle rice and late rice was consistent. The interannual change in ET of middle rice and late rice was not significant before 2003, and ET of middle rice and late rice showed a downward trend from 2003 to 2009, while ET of middle rice and late rice showed an upward trend from 2009 to 2017.

As can be seen from Figure 7, the interannual variation trend of precipitation is opposite to the interannual variation trend of ET overall in rice growing season, showing a negative correlation which is especially obvious in the early rice growth season. This may be that the study area is located in the subtropical monsoon region, and the east Asian summer circulation is enhanced during the early rice growth period, resulting in increased precipitation. ET can be regarded as one of the most important factors for indicating climate change at the catchment scale [55], and our research supports this view. Meanwhile, the interannual variation of ET is closely related with many ecological factors. This study's focus on the analysis of the time characteristics of ET and cannot directly explain the relationship between the interannual variation in ET and precipitation. If ET values are combined with meteorological factors, topographical factors, and vegetation changes, we may arrive at a more accurate and comprehensive understanding of the relationship between the interannual variation in ET and precipitation in the study area.

#### *3.4. Spatial Distribution of ET during Rice Growth in Different Hydrological Years*

Table 3 shows the precipitation and ET values during the rice growing season in different hydrological years, and the selection of a typical hydrological year is explained in Section 2.3. Almost no difference was observed between the ET of early rice in wet years and dry years, which might be due to the abundant precipitation of different hydrological years in the growing season of early rice. The ET of middle rice in different hydrological years showed significant differences, with a value of 393 mm in wet years, 423 mm in normal years, and 520 mm in dry years. The ET of late rice in wet years and normal years showed little difference, and the values in dry years were 330 mm and 367 mm. Overall, ET showed limited differences in the early and late rice in different hydrological years, but great differences in middle rice.

**Table 3.** ET of early, middle, and late rice during typical years in the GFPIS.

\* ET in 2008 could not be estimated because no images were available during the growing season of early rice.

Figures 8 and 9 show the ET frequency distribution and spatial distribution of rice in different hydrological years. In the wet year (2014), the early rice ET in the GFPIS was mainly within the range of 200–300 mm, with an average of 263 mm. The ET of the first, third, and fourth main canals of the irrigated area was between 263 mm and 283 mm. The early rice ET in the second main canal and the middle of the main trunk canal was lower than the average early rice ET in the whole irrigated area, and the early rice ET in the east trunk canal and the main canal was generally approximately 243–263 mm. In the dry year (2007), the ET of early rice was mainly in the range of 225–300 mm, with an average of 267 mm. Among them, the early rice ET in the first, second, and third main canal irrigated areas was higher than the average early rice ET value in the whole irrigated area, and the early rice ET in the fourth main canal and the end of the third main canal irrigated area was lower than the average value. These results indicate that the mean values of early rice ET do not differ significantly in different hydrological years, although significant differences occur in the spatial distribution of early rice ET.

The middle rice ET values in the wet year (2014) mainly ranged between 375 mm and 425 mm and the mean middle rice ET value in the whole irrigated area was 392 mm, with little spatial difference. In the normal year (2004), the middle rice ET mainly ranged from 360 mm to 450 mm and the average middle rice ET value in the whole irrigated area was 423 mm. In terms of spatial distribution, the middle rice ET in the east main trunk and the end of the fourth main canal was higher than 450 mm. In the dry year (2016), the middle rice ET value mainly ranged from 487–570 mm and the mean ET value was 525 mm in the whole irrigated area. There were differences in the spatial distribution of ET, with generally higher values at the end of the fourth main canal and the second main canal compared with the mean ET value in the whole irrigated area while the ET values in the middle section of the west main canal and the second main canal were lower than 525 mm.

**Figure 8.** Frequency distribution of rice ET in different hydrological years.

The late rice ET values in the wet year (2002) in the GFPIS mainly ranged from 300 mm to 350 mm, with an average value of 328 mm. In the western GFPIS, the late rice ET was generally higher than that of early rice and the ET was evenly distributed. The eastern GFPIS is a hilly area, and the late rice ET was generally lower than the mean value. In the normal year (2010), the late rice ET values in the GFPIS areas mainly ranged from 300 mm to 380 mm, with an average value of 347 mm. The late rice ET at the end of the second main canal and the fourth main canal was higher than that in the whole irrigated area (347 mm). The late rice ET in the third main canal and the first main canal irrigated area was higher than that in the whole irrigated area. In the dry year (2008), the late rice ET in the whole irrigated area mainly ranged from 325 mm to 400 mm, with an average value of 366 mm. The spatial distribution of the late rice ET was similar to that of the wet year, and the late rice ET in the western irrigated area was higher than the average value, while the ET in the eastern irrigated area was lower than the average value overall. Compared with the wet year, the ET in the whole irrigated area increased in dry years and mainly ranged from 40 mm to 60 mm.

**Figure 9.** Spatial distribution of rice ET in different hydrological years.

In general, the late rice ET in the western irrigated areas was higher than the mean late rice value in irrigated areas in different hydrological years. The late rice ET in the irrigated areas of the first and second main canals and the third canal was slightly higher than the mean late rice ET in the whole irrigated area in the three hydrological years. In the normal and dry years, the late rice ET at the end of the fourth main canal and the second main canal was much higher than the mean late rice ET value in the whole irrigated area. This finding is related to water resource dispatching measures in the GFPIS. In the case of drought, the first, second, and third main canals on the east main canal and west main canal should be adopted to temporarily control water inflow so that the water output can guarantee water for production and living in the middle and lower reaches of the irrigated area. The late rice ET of the eastern irrigated area and the hilly area near the west main water inlet in the three hydrological years was generally lower than the mean late rice ET value of the GFPIS [56]. Agricultural water consumption accounts for approximately 8% of the actual water consumption in the eastern irrigated area; therefore, the irrigation conditions in the irrigated area are not as good as those in the western irrigated area.

#### **4. Conclusions**

In this paper, the SEBAL model was applied to the subtropical climate region of southern China (GFPIS) to assess its rationality and applicability for daily scale estimates of paddy ET and to analyze the spatial and temporal distribution characteristics of early, middle, and late rice evapotranspiration in the GFPIS in different hydrologic years. The main conclusions are as follows:

Due to the large amount of water evaporation in paddy fields, the single-source SEBAL model needs further research to estimate ET in rice flooding periods, which will improve the accuracy of the SEBAL model in estimating ET in rice growing areas.

**Author Contributions:** Conceptualization, Y.L., J.C. and G.W.; Methodology, Y.L., J.C. and G.W.; Software, J.C. and G.W.; Validation, H.X (Hua Xie). and Y.Y.; Formal analysis, J.C. and G.W.; Investigation, J.C. and H.X. (Hengwang Xie); Resources, C.W. and Y.C.; Data curation, Y.Y. and C.W.; Writing—original draft preparation, J.C. and G.W.; Writing—review and editing, Y.L., G.W. and Y.C.; Visualization, G.W.; Supervision, Y.L.; Project administration, Y.L.; Funding acquisition, Y.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was financially supported by the NSFC-MWR-CTGC Joint Yangtze River Water Science Research Project (No. U2040213) and the National Natural Science Foundation of China (nos. 51979201 and 51779187).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Effects of Mulching on Maize Yield and Evapotranspiration in the Heihe River Basin, Northwest China**

**Qianxi Shen 1,2,3, Jun Niu 1,2,3,\*, Bellie Sivakumar 4 and Na Lu 1,2,3**

**Abstract:** Plastic film mulching is an effective way to manage agricultural fields in water shortage areas. Through increasing the soil surface temperature at the early stage of crop growth and reducing the soil evaporation during the whole growth period, plastic film mulching can realize the effect of water saving and yield increase. This study examined the effects of plastic film mulching on crop yield and evapotranspiration (*ET*) in the Heihe River basin in Northwest China. By using remote sensing data, the gridding G-AquaCrop model was built to simulate the maize yield and *ET* in the basin under conditions of film mulching and no-film mulching. Through an analysis of changes in maize yield and *ET* before and after film mulching, suitable areas for film mulching in the whole basin were identified. Through comparative analysis, it was found that after plastic film mulching, maize yield in 12–41% of the Heihe River basin increased to a certain extent, reaching 8%. Furthermore, film mulching decreased *ET* by 5–30% in 34–41% of areas planted with maize. Based on these results, suggestions were made on suitable areas for expansion of maize cultivation to balance the benefits of water saving and production increase with environmental pollution. Furthermore, the way of assessing the suitable mulching area is obtained by examining the meteorological condition directly. The results of this study are of great significance for rational allocation of agricultural production resources and efficient utilization of agricultural water resources.

**Keywords:** maize yield; evapotranspiration; plastic film mulching; G-AquaCrop model; Heihe River basin

#### **1. Introduction**

Plastic film mulching is a field management technique that results from advances in the field of material science. In addition to all the benefits of traditional biological mulching (i.e., to increase temperature, conserve soil moisture, and reduce evaporation) [1,2], plastic film mulching has several advantages, such as low cost, easy production, and good performance [3]. As a result, plastic film mulching has been widely used in agriculture in recent years [4–7]. It is likely that the use of plastic film mulching will continue to increase in the future, especially considering the increasing food demands due to population growth and the anticipated impacts of climate change.

In China, increasing population and climate change have already raised concerns about food security, and the situation is anticipated to only worsen [8–10]. Projections of more frequent and greater magnitude droughts [11–14] emphasize the need for water saving in the agricultural sector, especially since the sector is the largest consumer of water. In this regard, plastic film mulching offers an important means, as it has the benefits of water saving, especially with low cost and greater efficiency in production. Needless to

**Citation:** Shen, Q.; Niu, J.; Sivakumar, B.; Lu, N. Effects of Mulching on Maize Yield and Evapotranspiration in the Heihe River Basin, Northwest China. *Remote Sens.* **2022**, *14*, 700. https://doi.org/10.3390/rs14030700

Academic Editor: Pradeep Wagle

Received: 21 December 2021 Accepted: 31 January 2022 Published: 2 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

say, the benefits of plastic film mulching are largest when it is used in regions where water resources are already scarce and agriculture is a key sector, such as the Heihe River basin in Northwest China.

The Heihe River basin, located in the arid Hexi Corridor of Northwest China, is an important commodity grain production base in China [15,16]. Due to its geography and climatic conditions, this area has been facing a severe shortage of water resources for a long time [17]. Theoretically speaking, plastic film mulching has a great potential for application in this area. However, due to differences in climate and soil types over the region, not all areas in the basin can produce satisfactory results with film mulching [18]. Another issue is that the plastic film cannot be degraded in a short time after being discarded [19], which not only results in environmental pollution but also reduces crop yield to varying degrees. Jiang et al. (2017) used the dye tracer method to track soil water movement, and the results highlighted that the presence of residual plastic film fragments significantly influenced soil physical properties, altered soil water distribution, and decreased the matching degree of the flow distribution region and the maize root (densely rooted) zone [20]. Gao et al. (2019) used meta-analysis to show that the crop yield decreased significantly with increased time and residual plastic film exceeded 240 kg/ha [21]. Thus, the use of plastic film mulching without considering its applicability may aggregate negative effects. Therefore, in actual production, it is necessary to allocate agricultural production according to local conditions in order to realize efficient utilization of resources.

Producing more food with less water has been the goal of agricultural development in arid areas, especially where mulch is widely used. Crop water productivity (*WP*) refers to the ratio of crop yield to crop evapotranspiration (*ET*) during yield acquisition [22–24]. By definition, a larger *WP* is the goal pursued by researchers and farmers. For this reason, *WP* can be used as a key index to explore the potential and suitability of plastic film mulching in the Heihe River basin.

The AquaCrop model published by the Food and Agriculture Organization (FAO) is a typical water-driven crop growth simulation model [25–27]. It accurately simulates not only the crop yield but also the *ET* in the process of crop growth. Further, it divides the *ET* into soil evaporation (*E*) and crop transpiration (*Tr*). At the same time, compared with other crop growth simulation models (e.g., EPIC, DSSAT, WOFOST, APSIM, and WAVES) [28–32], this model requires fewer input parameters with stable performance and can simulate different field management practices, such as different irrigation methods, film mulching methods, and fertilization methods [25,26]. It is indeed a good choice to use the AquaCrop model to explore the effects of different field management practices on *WP* in arid areas. Ran et al. (2017) showed that the simulation performance of the AquaCrop model was better than that of the SIMDualKc model under deficit irrigation conditions through a 5-year field experiment data verification [33]. Compared with SIMDualKc model [34], the AquaCrop model not only considered the effect of soil water content, but also took into account the physiological damage (for example, a decrease in canopy coverage) caused by long-term water shortage when simulating *Tr*., Shen et al. (2019) used the AquaCrop model to simulate the changes in *ET* of seed maize with and without film mulching in the Shiyang River basin in Northwest China and verified the results with field data [35]. The results not only showed that the model had good simulation performance under both conditions but also indicated that plastic film mulching significantly reduced soil evaporation, thereby reducing water consumption during the growth period of the crop. The current AquaCrop model does not take into account the effects of plant diseases and pests on crop growth and yield. In this case, the simulation results are often the potential yield (defined by solar radiation and temperature) or attainable yield (limited by water and nutrient availability) rather than the actual yield (reduced by diseases, pests, and environmental stressors); this seems to be a common flaw in most crop growth models [36]. Tigkas et al. (2015) and Li et al. (2019) showed that when the independent variables were meteorological conditions or management measures, the crop model was an ideal tool. In this way, other influencing

factors can be set to have a neutral role in crop growth without affecting the interpretation of the results [37,38].

Despite these encouraging results, it is important to note that most simulation models such as AquaCrop are applicable only at the field scale. This presents certain limitations when applied to the Heihe River basin. This is because, due to the wide distribution of maize planting areas in the whole Heihe River basin, the spatial heterogeneities meteorological, soil and management variables in the region are large and are difficult to calculate and measure [39]. Therefore, although the AquaCrop model has been found to perform well for field-scale studies, it cannot be directly applied to regional studies due to limitations and difficulties in parameter estimation.

Advances in remote sensing technology have provided effective support for exploring the potential of plastic film mulching in the Heihe River basin using the AquaCrop model [40,41]. For example, many studies have proposed different important indexes in agricultural production using remote sensing data. Jin et al. (2020) and Toshihiro et al. (2006) estimated crop physiological indicators such as canopy coverage (*CC*) and growth period by using remote sensing data such as the enhanced vegetation index (*EVI*) and the three-band water index (*TBWI*) [42,43]. Lorite et al. (2013) and Han et al. (2019, 2020) used remote sensing technology to obtain *CC*, sowing date, aboveground biomass and other crop physiological indicators preliminarily. Then, the above indicators were introduced into the crop model and distributed simulation was carried out [44–46]. These results indicate that the coupling of remote sensing data and crop models can extend the application scale of existing crop models.

According to the above review, plastic film mulching has great application potential in arid areas of northwest China. However, few studies attempt to clarify the application effect and applicability of plastic film mulching in the whole Heihe River basin; this provided the motivation for the present study. In this study, the spatial distribution of water-saving and yield-increasing effects of plastic film mulching in the Heihe River basin was investigated and suggestions for suitable expansion areas were put forward. First, the maize sowing date and *CC* without statistical data were retrieved from *EVI* data and normalized difference vegetation index (*NDVI*) data, respectively. This lays a foundation for subsequent simulation and analysis. Due to the limited availability of land use data, only data for the period of 2011–2014 were considered. Further, combined with other products for the grids (i.e., meteorological data, land use data, soil variables), the gridding G-AquaCrop model was constructed. The G-AquaCrop model outputted maize yield and *ET* during the whole growth period under different scenarios (mulching/no-mulching). On the basis of reliable output results, the application effect and applicability of plastic film mulching of maize field in the Heihe River basin were quantified. At the same time, a new method to quickly obtain the applicability of plastic film mulching was proposed.

#### **2. Materials and Methods**

#### *2.1. Study Area*

This study focused on the maize cultivation area of the Heihe River basin (Figure 1) in the arid inland region of Northwest China. The basin is located in the middle of the Hexi Corridor, at 97◦42 ~102◦04 E, 37◦41 ~42◦42 N and covering an area of 142,900 square kilometers [47]. The maize cultivation areas in this study are mainly distributed in nine counties: Shandan, Minle, Ganzhou, Linze, Gaotai, Sunan, Suzhou, Jiayuguan, and Jinta in the middle and lower reaches of the Heihe River basin, as shown in Figure 1. In this area, the climate is mainly temperate continental arid climate; the annual average temperature is 4–8 ◦C; the annual average precipitation decreases from more than 250 mm in the east to less than 50 mm in the west; the annual average evaporation increases from the east to the west, from less than 2000 mm to more than 4000 mm; and the annual average sunshine duration reaches 3000–4000 h [48–50]. It is a relatively suitable area for agricultural development in the whole Heihe River basin. The main crops in the region are maize, spring wheat, spring barley, and spring rapeseed. Among these, maize occupies the largest area [51].

**Figure 1.** Location map of the Heihe River basin and the nine counties considered in this study.

*2.2. Acquisition of Spatial Heterogeneity Variables*

#### 2.2.1. Land Use Data

The land use data used in this study were from the Landsat remote sensing data obtained by the Heihe Plan Science Data Center (http://heihedata.org/ (accessed on 10 December 2019)), with a spatial resolution (*SR*) of 30 m. In order to matching other available forcing data, the 500 m *SR* is applied. The 500 m *SR* can also reflect spatial heterogeneity for the relatively flat agriculture land. The area cultivated under maize was re-samples to 500 m resolution using NEAREST by ArcGIS. This data was used as the basic data for this study and also as a reference for other data processing.

#### 2.2.2. Meteorological Data

The ERA (*ECMWF* re-analysis-interim) data of precipitation (P, mm), maximum air temperature (*Tmax*, ◦C), minimum air temperature (*Tmin*, ◦C), wind speed at 10 m height (*u*2, m/s), solar net radiation (*Rn*, MJ m−2 d−1), surface pressure (*SP*, kPa), and dew point temperature (*Tdew*, ◦C) was obtained from the European Centre for Medium-Range Weather Forecast (*ECMWF*, https://www.ecmwf.int/ (accessed on 11 August 2021)), with a spatial resolution of 0.125◦ (13.89 km) and a temporal resolution of one day. In this study, the spatial resolution was re-sampled to 500 m using NEAREST with the land use data as the reference.

The reference crop evapotranspiration (*ET*0, mm/d) was calculated by the Penman– Monteith method [52]:

$$ET\_0 = \frac{0.408 \cdot \Delta (R\_{\rm tr} - G) + \gamma \cdot \frac{900}{T + 273} \cdot \mu\_2 \cdot (\varepsilon\_s - \varepsilon\_a)}{\Delta + \gamma \cdot (1 + 0.34 \cdot \mu\_2)} \tag{1}$$

where *Rn* is the average net radiation of the crop surface; *G* is the soil heat flux, which is approximately equal to 0 in a day, *γ* is the hygrometer constant, *T* is the average temperature; *u*2 is the average wind speed at a height of 2 m, *es* is the saturated water pressure at average temperature, *ea* is the actual water pressure; and Δ is the slope of the water pressure curve.

#### 2.2.3. Soil Variables

The main soil types in the Heihe River basin include clay, sand, and silt. The distribution of soil types was obtained from the Digital soil mapping dataset of soil texture in the Heihe River basin (2012–2014) (spatial resolution was 100 m) [53] downloaded from the National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/zh-hans/ (accessed on 25 August 2020)). The distribution of the soil hydraulic properties (including soil saturated water content, field water capacity, and saturated water conductivity) were also derived from the Digital soil mapping dataset of hydrological features in the Heihe River basin (2012) (spatial resolution is 90 m) [54] downloaded from the same website. In this study, the spatial resolution was uniformly aggregated to 500 m by NEAREST.

#### 2.2.4. Sowing Date Retrieval

The sowing date is the starting point of time for crop growth and yield simulation by the AquaCrop model. Accurate acquisition of the spatial distribution data of sowing date in grids is of great importance for maize yield and *ET*. With the help of *MODIS-NDVI*, Han et al. (2019) obtained the sowing date of maize in the middle reaches of the irrigation area in the Heihe River through the threshold method [45]. However, many studies have shown that *EVI* has been improved to overcome the problems of atmospheric noise, soil background, and saturation on the basis of *NDVI*, and more and more studies choose to use *EVI* for phenological monitoring [55–59]. Therefore, in this study, MYDEV1D China 250 M *EVI* daily product of Geospatial Data Cloud (http://www.gscloud.cn/#page1/3 (accessed on 18 September 2020)) was used for sowing date retrieval. The sowing date was calculated using [60]:

> ⎧ ⎪⎪⎨ ⎪⎪⎩

$$\begin{array}{l}EVI\_{i} > \beta\\EVI\_{i+5} - EVI\_{i} > 0\\EVI\_{i+10} - EVI\_{i+5} > 0\\swing\ date = i - 10\end{array} \tag{2}$$

where *EV Ii* is the *EVI* value of day *i*, and *β* is the empirical threshold, which is closely related to crop variety and planting area. Due to the lack of empirical data, this study used the empirical sowing date data of the Yingke Irrigation District in the Heihe River basin, i.e., 22 April [61], to carry out the backward deduction so that the calculated sowing date generally falls on 22 April. Here, *β* is 0.22.

#### 2.2.5. Canopy Coverage Retrieval

Canopy coverage (*CC*) is an important variable connecting the maize plant body and the AquaCrop model, which reflects crop growth [62]. In this study, MYDND1M China 500 M *NDVI* monthly composite product of the Geospatial Data Cloud (http://www. gscloud.cn/#page1/3 (accessed on 23 September 2020)) was adopted. The spatial resolution of data was re-sampled to 500 m using NEAREST by ArcGIS. The CC was calculated using:

$$NDRVI = \frac{[(a+1)NDVI + (a-1)]}{[(a-1)NDVI + (a+1)]} \tag{3}$$

$$LAI = LAI\_{\max} \frac{WDRVI - WDRVI\_{\min}}{WDRVI\_{\max} - WDRVI\_{\min}} \tag{4}$$

$$\mathbb{CC} = 1.005 \times \left[1 - e \exp(-0.6LAI)\right]^{1.2} \tag{5}$$

where *WDRVI* is the wide dynamic threshold vegetation index [63], *LAImax* is the maximum leaf area index of maize, which was set as 5 [64]; and *α* is the empirical coefficient, which was set as 0.2.

#### *2.3. G-AquaCrop Model Parameterization & Yield and ET Simulation*

#### 2.3.1. G-AquaCrop Model

Since the AquaCrop model came into being, it has been widely used and recognized, especially in arid areas [25,26,33,35].

In the simulation of the AquaCrop model, the key process of conversion from *CC* to yield (*Y*) may be summarized as [26]:

$$\text{CC}^\* = 1.72 \text{CC} - \text{CC}^2 + 0.30 \text{CC}^3 \tag{6}$$

$$T\_I = K\_\mathcal{S} \mathcal{C} \mathcal{C}^\* K c\_{Tr,x} E T\_0 \tag{7}$$

$$B = K s\_b W P^\* \sum \frac{T\_r}{ET\_0} \tag{8}$$

$$Y = f\_{HI} H I\_0 B \tag{9}$$

where *CC*∗ is the adjusted canopy coverage, *KcTr*,*x* is the largest standard crop transpiration coefficient, *Ks* is the soil water stress coefficient, *Tr* is the plant transpiration, *WP*∗ is the normalized water productivity, *Ksb* is the air temperature stress coefficient, *H I*0 is the reference harvest index, and *fH I* is the water stress adjustment factor. In this study, water stress was not considered, and so both *Ks* and *fH I* were 1. According to Equations (5)–(8), accurate simulation of *CC* is a necessary prerequisite for the accurate simulation of maize *Y* and *ET* (For more detailed model mechanism process, please refer to Raes et al. (2009)).

However, since the Heihe River basin covers a large area and the meteorological, soil, vegetation, and other variables have significant spatial heterogeneity, it is impossible to use the point-scale AquaCrop model for simulation. Therefore, this study combined the advantages of AquaCrop plug-in v6.0, ArcMap 10.0, and Matlab R2018a software to build a G-AquaCrop model. The G-AquaCrop model uses ArcGIS 10.0 software to grid all spatial heterogeneity variables obtained by remote sensing. Each grid in the Heihe River basin was taken as an independent simulation unit and was simulated by the AquaCrop model. Then, the ArcGIS 10.0 software was used to present the simulation results of each grid in the basin map. Due to the large number of grids in the Heihe River basin, this study prepared each input file required for the model operation in advance and used Matlab R2018a to call AquaCrop plug-in v6.0 to simulate. Calendar Day was used for simulation in this study.

#### 2.3.2. Setting Simulation Scenarios for the G-AquaCrop Model

In order to compare the effects of plastic film mulching on *Y* and *ET* in the Heihe River basin, two simulation scenarios were set up in this study: (1) plastic film mulching; and (2) no plastic film mulching. In order to control the variables and highlight the influence of plastic film mulching, the irrigation schedule of the two scenarios was unified. Due to the lack of actual irrigation statistics, the irrigation schedule was set based on the irrigation experience of Yingke irrigation area. Yingke irrigation area is located in the middle reaches of the Heihe River basin, which could be regarded as a good representative for the region. The irrigation amount in Yingke irrigation area was 507–819 mm, and its average value (663 mm) was taken as the total amount of irrigation here. The irrigation events were four times, with border irrigation. The water stress was excluded from all scenarios in this study. It is assumed that there is no fertility stress, i.e., fertility is not limited. At the same time, due to the widespread application of pesticides in agricultural production for the study period, plant diseases and insect pests were not included in the two scenarios. In addition, according to the FAO recommended value, under the film mulching condition, the mulching proportion parameter was set to 70% and the mulching effect parameter was set to 100%. Under the condition of no film mulching, the film mulching proportion and the film mulching effect parameters were set as 0 (Table 1). Due to the lack of spatial distribution data of plastic film mulching in the Heihe River basin, this study considered the non-plastic film mulching scenario as the current basic scenario, and the simulation results of Y under the non-plastic film mulching condition were compared with the data in the 2012–2015

Statistical Yearbook of the Gansu Provincial Bureau of Statistics (http://tjj.gansu.gov.cn/ (accessed on 13 December 2020)). Table 1 shows the management parameter settings in different scenarios.

**Table 1.** Parameters for characterizing on-farm water management interventions in AquaCrop.

√: Management measures used in this study.

#### 2.3.3. Tuning and Simulation of G-AquaCrop Model

Increasing the accuracy of the simulation parameter in each grid can improve the simulation of *Y* in the G-AquaCrop model and the accuracy of *ET* simulation. In addition to the spatial heterogeneity variables that can be obtained from remote sensing data, there are some mechanical parameters involved in the G-AquaCrop model. According to past studies [35], the nine parameters in Table 2 are the most sensitive parameters about maize *Y* and *ET* in the model, which would be adjusted in this study.

**Table 2.** Optimization parameters and their default values and calibrated values.

As mentioned above, accurate simulation of *CC* is the premise of accurate simulation of maize *Y* and *ET*. Therefore, monthly *CC* values of remote sensing retrieval were taken as the verification data of each grid's parameter in the basic scenario (that is, without film mulching). With the help of Matlab R2018a software, a particle swarm optimization (*PSO*) algorithm was used to conduct gridding optimization and parameter tuning [65,66]. In order to highlight the effect of film mulching and simplify the simulation process, the same set of parameters was used for the corresponding grid of film mulching and nofilm mulching scenarios, and only the proportion parameter and effect parameter of film mulching were changed.

The objective of parameter optimization is to minimize the difference between the *CC* values obtained by remote sensing and which output by the model on the corresponding date. The objective function is:

$$\min f\_{\text{CC}} = \sqrt{\frac{1}{n} \sum (\text{CC}\_{RS} - \text{CC}\_S)^2} \tag{10}$$

where *CCRS* refers to the *CC* value of remote sensing retrieval, *CCS* refers to the *CC* value of the corresponding date output by AquaCrop model, and *n* is the number of *CCRS* during the growth period of maize in each grid.

Parameters were optimized via the server. The CPU was Intel (R) Xeon (R) Silver 4214 @2.20 GHz 2.19 GHz, and the operating memory was 128 GB. In the *PSO* algorithm, the number of particles and iterations were both set to 30. The elapsed time of optimization was 210 s for each grid. Due to the large number of grids in the research area, taking the year 2014 as an example, the number of grids was 9901. Therefore, all grids of each year were first grouped according to the sowing date to improve optimization efficiency. Some grids were selected for parameter optimization in each group, which (about 60) with the nearest grid number shared a set of parameter optimization results. The spatial distribution of parameter optimization objective function is shown in Figure 2, where, min*fCC* < 0.2 means the optimization result is Very good, 0.2 ≤ min*fCC* < 0.6 means the optimization result is Good, and min*fCC* > 0.6 means the optimization result is Normal. It can be seen from Figure 2 that satisfactory results were obtained during parameter optimization for 4 years, providing a solid foundation for the accurate simulation of *Y* and *ET*. When min*fCC* is reached, the input parameters of this model are adopted as the optimal parameters, as shown in Table 2.

After the parameters were adjusted, each grid-specific parameter was brought into the G-AquaCrop model to carry out the simulation of *Y* and *ET*. Through extensive literature reading, we obtained some ground experimental observations from published papers. The data were collected from three ground stations in the Heihe River basin, including maize yield and leaf area index (*LAI*). Using Equation (5), *LAI* was converted into *CC*. The data details are shown in Table 3. Then, we fitted the simulation results of the corresponding grid with the observation data to evaluate the simulation effect of G-AquaCrop model. Due to the limited amount of ground observation data, this study also averaged the *Y* per unit area of grid in the same county without film mulching and compared it with the *Y* per unit area of the county obtained through the Statistical Yearbook. Here, we used the following statistical evaluation indicators: the root mean squared error (*RMSE*), normalized root mean squared error (*NRMSE*), Willmott's index of agreement (*d*), and mean bias error (*MBE*).

$$RMSE = \left[\frac{\sum\_{i=1}^{n} (M\_i - S\_i)^2}{n}\right]^{0.5} \tag{11}$$

$$NRMSE = \frac{100}{\overline{M}} \left[ \frac{\sum\_{i=1}^{n} (M\_i - S\_i)^2}{n} \right]^{0.5} \tag{12}$$

$$d = 1 - \frac{\sum\_{i=1}^{n} \left(M\_i - S\_i\right)^2}{\sum\_{i=1}^{n} \left(|S\_i - \overline{M}| + |M\_i - \overline{M}|\right)^2} \tag{13}$$

$$MBE = \frac{1}{n} \sum\_{i=1}^{n} (S\_i - M\_i) \tag{14}$$

where *Mi* is field observation data or yearbook data, and *Si* is the simulated data. When *Mi* represents field observation data, n represents the quantity of observation data; when *Mi* represents yearbook data, *n* is the number of counties. For the *RMSE*, when the *RMSE* value is close to 0, the simulation results are in good agreement with the measured results. While for the *d*, the closer it is to 1, the better the simulation result is. The *NRMSE* can indicate the degree of overall consistency, and the *MBE* can represent positive and negative attributes of model bias.

**Figure 2.** Spatial distribution of parameter optimization objective function (min*fCC*) in the Heihe River basin for four different years: (**a**) 2011, (**b**) 2012, (**c**) 2013, and (**d**) 2014.

**Table 3.** Field observation data details.

#### *2.4. Relationship between Mulching Areas and Meteorological Variables*

The suitable areas for mulching were given by comparing the model output results under different conditions. Based on this result, the study attempts to explore the relationship between the mulching areas and their meteorological variables (accumulated temperature and precipitation). Here, the relationship between mulching areas and meteorological variables was expressed by a ratio, which of the number of grids recommended for mulching to the total number of grids under each meteorological condition:

$$r\_i = \frac{N\_{mulch}}{N\_{mulch} + N\_{no-mulch}} \tag{15}$$

where *i* refers to the specific total accumulated temperature value or total precipitation value in each grid during the whole crop growth period. When *i* stands for accumulated temperature, 1200 d°C < *i* < 3000 d°C; when *i* represents precipitation, 0 mm < *i* < 1000 mm. *Nmulch* is the number of grids mulching under condition of *i* meteorological value; *Nno*−*mulch* is the number of grids without mulching under *i* condition; *ri* is the ratio of the recommended grid number for film mulching to the total grid number under the same accumulated temperature value or precipitation value.

#### **3. Results**

#### *3.1. Spatial Distribution of Sowing Date*

Through the inversion method of the sowing date introduced above, the spatial distribution map of the maize sowing dates from 2011 to 2014 in the Heihe River basin was obtained, as shown in Figure 3. As can be seen from Figure 3, the sowing date of maize in the Heihe River basin was approximately from the 91st day to the 142nd day of a year, with a span of 52 days. In most areas, the sowing dates are concentrated on the 111th–125th day, around 22 April. From the spatial perspective, this part of the region is densely distributed in the Ganzhou and Linze counties. The maize areas in these two counties are irrigated areas, where unified field management measures are implemented. The sowing date of maize in the irrigated areas is about 22 April, which is consistent with the results obtained in this study. However, in areas outside the irrigated areas, such as Jinta and Suzhou, there are certain differences between the sowing dates and the irrigated areas due to differences in geography, climate, and crop varieties. It is obvious that the sowing date in the northern part of the Ganzhou region is generally early, which may be related to the maize varieties planted in this region. Ganzhou is the largest maize seed production base at the county level in China. However, some studies have shown that in actual production, the seed maize should be separated from the sowing date of field maize with a flexible time of 30–40 days [69] to better carry out seed production. This is why the sowing date in the northern region of Ganzhou is regularly advanced.

**Figure 3.** Sowing date of maize in the Heihe River basin for four different years: (**a**) 2011, (**b**) 2012, (**c**) 2013, and (**d**) 2014.

#### *3.2. Spatial Distribution of Y under Different Scenarios*

Figure 4 shows the spatial distribution of maize *Y* per unit area in the Heihe River basin from 2011 to 2014 under the two scenarios. Figure 4a,c,e,g represents the no-film mulching scenario for the four respective years (top to bottom), and Figure 4b,d,f,h represents the film mulching scenario. In general, the maximum *Y* per unit area of maize in the Heihe River basin was 17.6–20.7 t/ha, while the average *Y* was only 8.6–10.0 t/ha. These values indicated that there was a large difference in the spatial distribution of maize *Y* per unit area. At the same time, after film mulching, the maximum *Y* did not increase significantly, but in a local range, film mulching did improve the *Y* of maize to a certain extent. From one point of view, plastic film mulching is already in use in some areas. If plastic film mulching was covered in a wide range, there would be two consequences: first, plastic film mulching could not reach the desired effect; second, it would put great pressure on the ecological environment.

**Figure 4.** Spatial distribution of simulated maize yield in the Heihe River basin from 2011 (top) to 2014 (bottom): (**a**,**c**,**e**,**g**) are the *Y* values without mulching; and (**b**,**d**,**f**,**h**) are the *Y* values with mulching.

Figure 5 shows the comparison between the observed and simulated *Y* and *CC* of maize at the three field stations in the Heihe River basin. Through consulting the published literature, we obtained some data of maize planting experiment in Heihe River basin. After excluding the data with inconsistent planting years and unclear planting locations, data of the three field stations in Table 3 were finally selected for comparison. The *d* and *R*2 values for maize *Y* and *CC* were larger than 0.90, indicating good agreement for both *Y* and *CC*. The simulated *Y* slightly underestimated the observed values, with *MBE* of −0.51 t/ha. While the simulated *CC* slightly overestimated the observed values, with *MBE* of 1.6%. Considering all the statistical indicators, the G-AquaCrop model has good simulation performance. The results of Zhang et al. (2021) also included 5 *Y* values [70]. All the statistical indicators about *Y* in our study are superior to the results of Zhang et al.

**Figure 5.** Comparison between the observed and simulated (**left**) *Y* and (**right**) *CC* of maize at the three field stations in the Heihe River basin.

Table 4 shows the fitting effect of maize *Y* per unit area at the county level in the Heihe River basin from 2011 to 2014. The results showed that, under the condition of no-film mulching, the *RMSE* and *NRMSE* of spatial average *Y* reached 1.29 t/ha and 14.37 (dimensionless), respectively. The absolute values of *MBE* were not larger than 0.77 t/ha. However, the absolute value of maize yield *MBE* was 1.19 t/ha in the study of Zhang et al. (2021) [70]. For 9901 grids in the Heihe River basin, the fitting accuracy of our study is acceptable. After film mulching, the *Y* in each county was generally improved, to a certain extent. However, the simulated *Y* in Shandan in 2013 and 2014 was abnormally low. According to the meteorological data of these two years, spring drought of different degrees occurred in Shandan from January to April in 2013 and 2014. The government and agricultural management departments issued early warning information in advance and suggested farmers to take countermeasures. Therefore, the Statistical Yearbook data on *Y* did not show much decline. However, in the G-AquaCrop model simulation, empirical irrigation data was still used, resulting in abnormally low values. Similarly, spatial averaging weakened the local effect of film mulching, as will be described in detail in Section 3.3.

**Table 4.** Maize yield simulation in counties in the Heihe River basin for four different years (2011–2014). (YY represents the yearbook yield per unit area, YN represents the yield per unit area of the model without film mulching, and YM represents the yield per unit area of the model without film mulching in t/ha.).

In general, the distribution of maize in the irrigated areas of Ganzhou and Linze was relatively concentrated and the *Y* value was relatively stable, which is closely related to the mature planting system and field management measures in the irrigated areas. Han et al. (2019) obtained the spatial distribution of maize *Y* in the Yingke Irrigated Area in 2012 and 2015 [45]. In terms of quantity, the maize *Y* in the Yingke Irrigated Area ranged from 7.8 to 13.7 t/ha, which was consistent with the results obtained in this study. At the same time, Han et al. (2019) also found that the reason why the *Y* of 2015 was higher than that of 2012 might be that the sowing date of 2015 was about 10 days earlier than that of 2012. Compared with Figures 3 and 4, in this study, the sowing date of northern Suzhou was earlier and its *Y* was higher, which was consistent with the results of Han et al. (2019).

#### *3.3. Change in Y after Film Mulching*

The Figure 6 shows the spatial distribution of *Y* change ratio after plastic film mulching in the Heihe River basin from 2011 to 2014. It can be clearly seen from the spatial threedimensional diagrams that the *Y* change of most areas in the basin is −2–0%, and the negative value mainly comes from the random error of the model. In this study, the maize *Y* in the irrigated area after film mulching did not increase significantly, while the *Y* in Yumen, Jiayuguan, Jinta, Suzhou, Sunan, Shandan, and Minle counties significantly increased after film mulching, with the maximum increase of nearly 8%. There are two main reasons. First, the *Y* of maize in the irrigated area hardly increased after plastic film mulching, which may be due to the improvement of irrigation system, water and fertilizer management, and field management measures in the irrigated area at the present stage. Plastic film mulching may have been used in a large range, and the *Y* potential of maize has been fully tapped. However, since the spatial distribution of the film-mulched areas in the Heihe River basin could not be obtained at the present stage, this study assumed that the whole Heihe River basin was not film-mulched during parameter adjustment, and the *CC* data of remote sensing retrieval used in parameter adjustment was the actual *CC*. It is possible that the *CC* used for the parameter adjustment may be relatively larger in the irrigated area. Second, Jiayuguan, Jinta, and Suzhou are located in the middle and lower reaches of the Heihe River basin, whose latitudes are higher than that of the other counties, so the corresponding temperature is relatively low and the annual average precipitation is also less. Film mulching can increase soil surface temperature and reduce soil evaporation, thus promoting various physiological processes of maize plants. Therefore, after film mulching, the maize *Y* of Jiayuguan, Jinta, and Suzhou improved to a certain extent.

**Figure 6.** Spatial distribution of maize *Y* improvement degree after mulching and the proportion of each improvement degree in the Heihe River basin for four different years: (**a**) 2011, (**b**) 2012, (**c**) 2013, and (**d**) 2014.

#### *3.4. ET Changes after Film Mulching*

Figure 7 shows the spatial distribution of changes in *ET* after plastic film mulching in the Heihe River basin from 2011 to 2014. As can be seen from these diagrams, after plastic film mulching, *ET* in the whole maize planting area in the Heihe River basin decreased, with the maximum reduction reaching 30%. Further, *ET* in most areas decreased between 0.9 and 10.0%, including the irrigated areas with relatively perfect field management measures. The areas with large *ET* reduction were concentrated in Shandan and Minle, especially in Shandan, where the *ET* reduction reached 20–30% in 2014. This was mainly due to a large degree of spring drought in Shandan from January to April 2014, which reduced the soil evaporation by plastic film mulching.

**Figure 7.** Spatial distribution of maize *ET* reduction degree after mulching and the proportion of each reduction degree in the Heihe River basin for four different years: (**a**) 2011, (**b**) 2012, (**c**) 2013, and (**d**) 2014.

#### *3.5. Suggestions for Film Mulching According to Changes in Yield*

On the basis of the spatial distribution in Figure 6, a statistical chart of *Y* change range after plastic film mulching in the Heihe River basin from 2011 to 2014 is shown in the bar chart of each plot in Figure 6.

As can be seen, in four years, the areas where *Y* barely changed after plastic film mulching accounted for 59.8–87.8% of the maize cultivation area in the entire Heihe River basin. Among the areas where *Y* increased, the areas where *Y* increased by 0.1–5.0% accounted for the largest proportion, and the areas where *Y* increased by 5.1–8.0% accounted for less than 0.9% of the maize planting area. As mentioned above, the areas with no significant change in maize *Y* after plastic film mulching was mainly the irrigated areas in the middle reaches of the river basin, where various field management measures are relatively mature and crop *Y* potential is largely developed. On this basis, the areas with improved *Y* in the simulation should be encouraged to expand the area mulched with film in the actual production at the present stage to increase the maize yield. Figure 8 presents the proposed map of the mulched area according to different targets, including changes of maize *Y* (bottom layer in each plot), *ET* (middle layer), and *Y* and *ET* (top layer). As can be seen from the bottom layer in each plots, it is recommended to continue to expand the area of plastic film mulching in the southern areas of Jinta, Suzhou, Minle, and Shandan to obtain a higher yield. Irrigated areas continue to maintain the current status of field management.

**Figure 8.** The recommendations for mulching according the change of maize *Y* (bottom map in each figure); *ET* (middle map in each figure); and maize *Y* and *ET* (top map in each figure) after mulching for four different years: (**a**) 2011, (**b**) 2012, (**c**) 2013, and (**d**) 2014. The red area is recommended for mulching, and the yellow area is not recommended for mulching.

#### *3.6. Suggestions Based on ET Changes*

On the basis of the spatial distribution in Figure 7, we summarized the statistical chart of changes in *ET* after plastic film mulching in the Heihe River basin from 2011 to 2014, as shown in the pie chart of Figure 7. From 2011 to 2014, after plastic film mulching, *ET* in the whole Heihe River basin decreased by varying degrees, while 59.8–66.6% of the region only decreased by 5.0% (the pie chart of Figure 7). Considering that the plastic film cannot be degraded in a short time after use, leading to white pollution to the ecological environment of the area, the water-saving benefits brought by the plastic film mulching will be greatly reduced. Therefore, it is not necessary to expand film mulching in areas where *ET* reduction is less than 5%. Increasing the plastic film mulching area in the areas where *ET* reduction exceeds 5% can reduce the water consumption of crops to a certain extent, thus promoting the sustainable development of water-saving agriculture. With the above results, a proposed film mulching map for the Heihe River basin was drawn based on the distribution of *ET* reduction during 2011–2014, as shown in the middle layer of each plot in Figure 8. The results indicate that it is more suitable to increase the area of film mulching in Ganzhou, Minle, Shandan, and southeast of Linze in the Heihe River basin, while it is not necessary to increase it in the other areas.

#### *3.7. Suggestions Based on Y and ET Changes*

The ability of plastic mulch to reduce soil evaporation and increase yield provides a feasible solution to improve crop water productivity. Therefore, combined with the lower and middle layers of each plot in Figure 8, it is recommended to apply film mulching in areas with increased *Y* and reduction of more than 5% in *ET*, mainly distributed in some areas of Ganzhou, Minle, Shandan, and other counties. Mulching film in this area can improve *WP* and achieve water saving and yield increase (the upper layer of each plot in Figure 8).

#### *3.8. Relationship between Mulching Areas and Meteorological Variables*

Relationship between mulching areas and meteorological variables were studied during the period 2011–2014 (Figure 9). As can be seen from Figure 9a, the proportion of film mulching is the highest when the accumulated temperature is between 2100 and 2400 d◦C. When the accumulated temperature exceeds 2700 d◦C and below 1800 d◦C, there is almost no film mulching. The reason is obvious: due to full irrigation, a large accumulated temperature can support the growth of maize independently, while too small accumulated temperature cannot meet the physiological requirements of maize germination. Even with measures such as film mulching, maize cannot grow normally. Figure 9b shows the relationship between film mulching areas and precipitation. As the study area is located in an arid area, the highlights in the figure are concentrated in the part with low precipitation. It can be seen that the proportion of film mulching is highest when the precipitation is between 300 and 400 mm. Figure 9c shows the comprehensive influence of accumulated temperature and precipitation. With the increase of accumulated temperature and the decrease of precipitation, the proportion of film mulching increases. Meteorological conditions suitable for mulching are concentrated in the lower right corner of Figure 9c. By referring to this figure, you can clearly know the accumulated temperature and precipitation most suitable for plastic film mulching. Conversely, when you want to know whether the climatic conditions of an area are suitable for mulching, you can also easily obtain the results by referring to it. Therefore, Figure 9c provides a new method for judging whether mulching is suitable for each maize cultivation area in the Heihe River basin, which greatly improves the efficiency of determining.

**Figure 9.** *Cont*.

**Figure 9.** Proportion of mulching under different meteorological conditions in the Heihe River basin for four years: (**a**) accumulated temperature, (**b**) precipitation, (**c**) accumulated temperature and precipitation.

#### **4. Discussion and Conclusions**

To study the effects of plastic film mulching on maize yield and *ET* in arid areas, we chose the Heihe River basin as a typical representative to explore the impacts at the regional scale. The gridding large-scale G-AquaCrop model was established fully considering spatial heterogeneity by using remote sensing data. The parameters of the G-AquaCrop model were optimized by the particle swarm optimization method. The yield and canopy coverage of maize were verified by field data, and the *R*2 was larger than 0.92.

There is a large amount of plastic-mulched farmland in China [71]. However, the records of the current situation of plastic film mulch used on all statistical yearbooks only included the mulching weight used. The spatial distribution of plastic film mulch was not recorded. This would undoubtedly cause some obstacles to exploring the application effect and applicability of plastic film mulching in the region. Some scholars used remote sensing technology to identify the situation of plastic film mulch [72–74]. They often obtained original images from UVA or high-resolution satellites. For a large area, a large amount of input is required. At the same time, the crop canopy coverage was extensive when it reached the middle and late growth period. The crop canopy would cover the ground. This could cause the identification results to be not reasonable.

The skillful strategy in our modeling approach realizes relatively reliable simulation in evaluating the effects of film mulching in the studied region. We set up two scenarios: no-mulching and plastic film mulching, and the comparison of the difference between the two scenarios is focused. There were two possible situations if there was little difference between the two scenarios. In the first situation, the area was already mulched with plastic film. The second situation was that plastic film mulch was not applicable in the area. In either situation, the area did not need additional mulching. Too much plastic film mulching not only alters soil moisture distribution, but also reduces crop yields over time. As mentioned earlier, the G-AquaCrop model does not take into account crop disease and insect infestation. Meanwhile, as pesticides have been widely used in the Heihe River basin, the effects of disease and pests are not considered here. From the perspective of the model mechanism, the limitation caused by field data loss was avoided by comparing the two scenarios.

Compared with no mulching, plastic film mulching could retain soil moisture and significantly reduce soil evaporation. The biggest drop is 30%. Plastic film mulching could improve maize yield, but the effect of it on maize yield varied spatially. Considering the environmental pressure caused by plastic film mulching, the suggestions for mulching expansion were given from different angles. In general, expanding the plastic film mulching area in the upper and lower reaches of the Heihe River basin could produce greater positive benefits. The middle reaches of the Heihe River basin where the irrigation area is located have been fully explored. This research provides a quantitative and spatial assessment of plastic film mulching application, which helps better to understand management practices.

There are also some limitations of the current study work. Although the results of this study by comparing the differences between the two scenarios are reliable, they can be further improved. With the emergence of additional higher-resolution remote sensing products, we will continue to carry out remote sensing retrieval of actual film mulching conditions and bring data into our research to obtain more accurate and reliable simulation results. The study period is relatively short due to the available land use data on spatial crop varieties, which is important for the regional simulation study. The used data for the period of 2011–2014 was well evaluated by using Google Earth high-resolution images and field survey data, and the overall accuracy reached 92.19%. Meanwhile, the 4 years for the period of 2011–2014 vary in terms of the precipitation variability. For example, the year 2013 and 2014 is the spring drought year, and the year 2012 is with abundant precipitation. The regional precipitation variability has been reflected to a certain degree. Long-term data support would be expected for the next-step investigation.

In terms of research results, the G-AquaCrop model has the original advantages of the AquaCrop model, while overcoming the limitation that the AquaCrop model can only be used on field points. Due to the water-driven attribute of the G-AquaCrop model, its outputs were sensitive to water. Ganzhou, Linze, Gaotai, Sunan, and Suzhou counties are located in the middle reaches of Heihe River basin, which are closed to the Yingke irrigation area. The irrigation schedule in this study was set according to the experience of Yingke irrigation area. Therefore, the yield simulation results in these counties were relatively good. As can be seen from Table 4, the G-AquaCrop model only slightly overestimated the maize yield of these counties. This was also a reflection of the full irrigation (no water stress) setup in this study. Jiayuguan, Jinta, and Yumen counties are located in the lower reaches of the Heihe River basin and have less precipitation. The maize yields of the three counties in Table 4 were underestimated. However, it also reflected that the G-AquaCrop model could well capture water changes, which was suitable for the northwest arid region with water shortage.

Moreover, based on mathematical statistics, the accumulated temperature and precipitation conditions with the maximum probability of plastic film mulching were obtained. The result provided a new method for quickly judging whether a certain area was suitable for plastic film mulching.

**Author Contributions:** Conceptualization, J.N.; methodology, J.N. and Q.S.; software, Q.S.; validation, N.L.; formal analysis, Q.S.; investigation, Q.S.; resources, J.N.; data curation, Q.S. and N.L.; writing—original draft preparation, Q.S.; writing—review and editing, J.N., N.L. and B.S.; visualization, Q.S.; supervision, J.N.; funding acquisition, J.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was financially supported by the National Natural Science Foundation of China (51621061, 51679233), the Program of Introducing Talents of Discipline to Universities (B14002) provided by Ministry of Education of the People's Republic of China.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data used in this research will be available (by the corresponding author), upon reasonable request.

**Acknowledgments:** Thanks for land use data provided by the Heihe Plan Science Data Center; meteorological data provided by the European Centre for Medium-Range Weather Forecast; soil types and hydraulic properties data provided by the National Tibetan Plateau Data Center; MODIS-EVI/NDVI data provided by Geospatial Data Cloud. The authors thank the reviewers for their constructive comments and useful suggestions on earlier versions of this manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Multiple Effects of Topographic Factors on Spatio-Temporal Variations of Vegetation Patterns in the Three Parallel Rivers Region, Southeast Qinghai-Tibet Plateau**

**Chunya Wang 1,2,3, Jinniu Wang 1,2,\*, Niyati Naudiyal 1, Ning Wu 1, Xia Cui 4, Yanqiang Wei 5 and Qingtao Chen 3**

**Abstract:** Topographic factors are critical for influencing vegetation distribution patterns, and studying the interactions between them can enhance our understanding of future vegetation dynamics. We used the Moderate-resolution Imaging Spectroradiometer Normalized Differential Vegetation Index (MODIS NDVI) image dataset (2000–2019), combined with the Digital Elevation Model (DEM), and vegetation type data for trend analysis, and explored NDVI variation and its relationship with topographic factors through an integrated geographically-weighted model in the Three Parallel Rivers Region (TPRR) of southeastern Qinghai-Tibet Plateau (QTP) in the past 20 years. Our results indicated that there was no significant increase of NDVI in the entire basin between 2000–2019, except for the Lancang River basin. In the year 2004, abrupt changes in NDVI were observed across the entire basin and each sub-basin. During 2000–2019, the mean NDVI value of the whole basin increased initially and then decreased with the increasing elevation. However, it changed marginally with variations in slope and aspect. We observed a distinct spatial heterogeneity in vegetation patterns with elevation, with higher NDVI in the southern regions NDVI than those in the north as a whole. Most of the vegetation cover was concentrated in the slope range of 8~35◦, with no significant difference in distribution except flat land. Furthermore, from 2000 to 2019, the vegetation cover in the TPRR showed an improving trend with the changes of various topographic factors, with the largest improvement area (36.10%) in the slightly improved category. The improved region was mainly distributed in the source area of the Jinsha River basin and the southern part of the whole basin. Geographically weighted regression (GWR) analysis showed that elevation was negatively correlated with NDVI trends in most areas, especially in the middle reaches of Nujiang River basin and Jinsha River basin, where the influence of slope and aspect on NDVI change was considerably much smaller than elevation. Our results confirmed the importance of topographic factors on vegetation growth processes and have implications for understanding the sustainable development of mountain ecosystems.

**Keywords:** MODIS NDVI; vegetation cover; trend analysis; Sen + Mann-Kendall; topography; GWR; Qinghai-Tibet Plateau

#### **1. Introduction**

Being an important component of terrestrial ecosystems, vegetation plays a crucial role in maintaining ecosystem stability, functioning, and services [1,2]. The dynamics of vegetation patterns is a complex and prolonged process, driven by multiple biotic, abiotic, and

**Citation:** Wang, C.; Wang, J.; Naudiyal, N.; Wu, N.; Cui, X.; Wei, Y.; Chen, Q. Multiple Effects of Topographic Factors on Spatio-Temporal Variations of Vegetation Patterns in the Three Parallel Rivers Region, Southeast Qinghai-Tibet Plateau. *Remote Sens.* **2022**, *14*, 151. https://doi.org/ 10.3390/rs14010151

Academic Editors: Jingzhe Wang, Zhongwen Hu, Yangyi Wu and Jie Zhang

Received: 7 December 2021 Accepted: 26 December 2021 Published: 30 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

anthropogenic factors such as climate change, land use change, and ecological engineering measures amongst others [3]. Being one of the most fragile terrestrial ecosystems, mountainous ecosystems are rather sensitive to global climatic change and human activities. A large vertical gradient and unique topographical conditions make mountainous regions the most abundant land unit on the earth and a key area for global biodiversity conservation [4,5]. Vegetation cover change of mountainous ecosystems have a much more significant and sensitive feedback on climate and other factors than plain areas. Therefore, the analysis of the process and driving factors of mountain vegetation dynamics has always been a key concern amongst global climatic change and mountain research community [6,7].