Contribution of Biophysical Factors to Regional Variations of Evapotranspiration and Seasonal Cooling Effects in Paddy Rice in South Korea

: Previous studies have observed seasonal cooling effects in paddy rice as compared to temperate forest through enhanced evapotranspiration ( ET ) in Northeast Asia, while rare studies have revealed biophysical factors responsible for spatial variations of ET and its cooling effects. In this study, we adopted a data fusion method that integrated MODIS 8-day surface reﬂectance products, gridded daily climate data of ground surface, and a remote sensing pixel-based Penman-Monteith ET model (i.e., the RS–PM model) to quantify ET patterns of paddy rice in South Korea from 2011 to 2014. Results indicated that the regional variations of the rice-growing season ET ( RGS-ET , the sum of daily ET from the season onset of rapid canopy expansion ( SoS ) to the end of the rice-growing season ( EGS )) were primarily inﬂuenced by phenological factors (i.e., the length of growing period-LGP ), followed by growing season mean climatic factors (i.e., vapor pressure deﬁcit- VPD , and air temperature). For regional variations of the paddy ﬁeld ET ( PF-ET , the sum of daily ET from the ﬁeld ﬂooding and transplanting date detected by satellite observations ( FFTD sat ) to SoS , and to EGS ), the extents were substantially reduced, only accounting for 54% of the RGS-ET variations. The FFTD sat and SoS were considered critical for the reduced PF-ET variations. In comparison to the temperate forest, changes in monthly ground surface air temperature ( T s ) in paddy ﬁelds showed the V-shaped seasonal pattern with signiﬁcant cooling effects found in late spring and early summer, primarily due to a large decline in daytime T s that exceeded the nighttime warming. Bringing FFTD sat towards late spring and early summer was identiﬁed as vital ﬁeld management practices, causing signiﬁcant declines in daytime T s due to enhanced ET . Results highlighted climate-warming mitigation by paddy ﬁelds due to early ﬂooding practices. estimates (mm m − 2 day − 1 ). The background environment in scenario 1 included incident solar radiation ( R s ↓ ) of 19.0 MJ m − 2 day − 1 , daily relative humidity ( RH ) of 70% and the ﬂooded bare land. Other model parameters were set at optimal levels except for the one used for the speciﬁc scenario analysis. In scenario 2, to evaluate the effects of canopy resistance ( r s , s m − 1 ), R s ↓ was set at 19.0 MJ m − 2 day − 1 , with daily RH of 70%, surface albedo of 0.18, and leaf area index ( LAI ) of 4.0 m 2 m − 2 . In scenario 3, there were the same climate variables as in scenario 1 but with variying K cb values.


Introduction
The FAO's annual proceedings have reported that around 167 million hectares of global land were plowed for paddy rice plantation in 2017, and more than 90% of global paddy rice areas were distributed in Northeast and Southeast Asia [1]. Rice is conventionally grown in fields that are flooded at different depths during the land preparation period (i.e., the first date of field flooding and transplanting detected by satellite observations-FFTD sat , herein referred to as the FFTD sat nominated surrogate of paddy rice land preparation, as commonly seen in the remote sensing research on regional paddy rice) until the end of the rice-growing season (EGS). Paddy rice with a water layer on the soil surface develops a wetland environment, which increases the water availability for evapotranspiration (ET). As a major component of surface energy balance, ET can facilitate the removal (1) Field-to-field changes in the FFTD could become an important biophysical factor influencing spatial and temporal variations of ET among paddy fields. (2) Considering the temperate forest being adjacent to paddy rice planting areas, the monthly T s of paddy fields may not always be higher than that of the temperate forest.

Materials and Methods
The pathway diagram of earth observations and numerical simulations to test the two hypotheses through analyzing the correlations among geographical variations of the ET, T s , and biophysical factors of paddy field is illustrated in Figure 1. The flooding single classification method using the Moderate Resolution Imaging Spectroradiometer 8-day surface reflectance products (MOD09A1) was adopted to detect the extent of paddy rice fields and quantify the FFTD of paddy field in South Korea from 2011 to 2014. Potential paddy rice pixels, detected by the flooding single classification approach [29], were geographically overlapped by a 5m resolution land use and land cover map (Supplementary Materials SP1 and Figure S1); such pixels with homogeneity percentage ≥90% (i.e., rice area of an individual grid cell was ≥90%) were filtered. Daily vegetation indices (VIs) per paddy pixel were produced using the double logistic model, which help develop daily LAI and K cb per paddy pixel according to field relationships among the VI of interest, LAI, and K cb . Geographical variations of the ET of paddy fields were obtained via fusing daily LAI, K cb , and meteorological variables to drive the RS-PM model that was comprehensively calibrated and validated using eddy covariance observations. Phenological metrics of paddy rice were extracted from the MODIS VIs data using the TIMESAT (i.e., a software package for analysing time-series of satellite sensor data) that was validated by field observations. The detailed processes involved by earth observations and numerical simulations are stated in the following sections.

Mapping the Extent of Paddy Rice Using MOD09A1 and National Census Datasets
The online repository of MOD09A1 product from 2011 to 2014 provides the information about the quality assurance (QA) for each geo-reference correcting surface reflectance imagery, indicating the presence of cloud cover and other variables. To minimize the error associated with cloud or aerosol contamination, per-pixel data with poor QA and/or blue reflectance of ≥0.2 were considered as outliers [29]. Such abnormal data were replaced with the mean of the subsequent and previous reflectance datasets. Reflectance data in blue (ρ blue , 459-479 nm), red (ρ red , 620-670 nm), near-infrared (ρ nir , 841-876 nm) and shortwave infrared (ρ swir , 1628-1652 nm) bands were used to calculate VIs, i.e., the enhanced vegetation index (EVI) [30], the normalized difference vegetation index (NDVI) [31], the land surface water index (LSWI) [29], and the optimized soil-adjusted vegetation index (OSAVI) [32]. EVI is the vegetation index sensitive to vegetation greenness, whereas LSWI is highly sensitive to land surface moisture content.

Mapping the Extent of Paddy Rice Using MOD09A1 and National Census Datasets
The online repository of MOD09A1 product from 2011 to 2014 provides the information about the quality assurance (QA) for each geo-reference correcting surface reflectance imagery, indicating the presence of cloud cover and other variables. To minimize the error associated with cloud or aerosol contamination, per-pixel data with poor QA and/or blue reflectance of ≥0.2 were considered as outliers [29]. Such abnormal data were replaced with the mean of the subsequent and previous reflectance datasets. Reflectance data in blue (ρblue, 459-479 nm), red (ρred, 620-670 nm), near-infrared (ρnir, 841-876 nm) and short-wave infrared (ρswir, 1628-1652 nm) bands were used to calculate VIs, i.e., the enhanced vegetation index (EVI) [30], the normalized difference vegetation index (NDVI) [31], the land surface water index (LSWI) [29], and the optimized soil-adjusted vegetation index (OSAVI) [32]. EVI is the vegetation index sensitive to vegetation greenness, whereas LSWI is highly sensitive to land surface moisture content. It is required to distinguish permanent water bodies from seasonally flooded rice fields. We defined a pixel covered by a water body when the year-round VIs value of that pixel had both the NDVI < 0.10 and LSWI [29]. Rice cropping systems are practiced primarily on great flat plains in the western part of South Korea. There are considerable changes in the elevation of the forested lands surrounding rice plains. A digital elevation map (GTOPO30, Global 30-Arc Second Elevation Dataset, https://rda.ucar.edu/datasets/ds758.0/, accessed date: 21 October 2020) was used to retrieve the slope of each pixel by the ArcGIS 10.2 (Environmental Systems Research Institute, Redlands, CA, USA). Areas with a slope of over 10% were considered as non-paddy rice zones.
In a previous study, Xiao et al. [29] proposed a flooding single classification method (LSWI + T ≥ EVI, T is a threshold parameter) with a constant T of 0.05 to identify paddy rice areas and evaluate the FFTD sat [33]. They designated the first date of the EVI 8-day period meeting the flooding single classification approach at a paddy pixel as the FFTD sat of the paddy pixel. Paddy pixels, obtained by the flooding single classification method, were denoted as potential paddy pixels, as they may still contain numerous small landholdings of subpixels. We assumed that the land cover and land use type (agriculture, forest, and grassland) throughout South Korea were relatively stable from 2011 to 2014 ( Figure S2). A spatial information system embedded in the ArcGIS 10.2 software was used to stack the Remote Sens. 2021, 13, 3992 5 of 22 5-m land cover map and the map of potential paddy pixels, aiming to count the number of paddy rice fields within each potential paddy pixel. The 5m spatial resolution land cover map of the raster form was developed by the Ministry of the Environment Republic of Korea through the comprehensive county field census and remote sensing retrieval by visual interpretation (For detailed information, please refer to Jeong et al., 2012 [18] and Figure S1). The Figure S3 in Supplementary Materials SP1 shows an example of the selection of homogenous paddy pixels in the Charyeong rice plain. Potential paddy pixels with homogeneity of ≥90% throughout South Korea were regarded as pure paddy pixels. Meanwhile, we assumed that the influences of other land-use types that were less than 10% in each pure paddy pixel would be randomly detrended. The same method for selecting pure paddy pixels was used to identify pure forest pixels. The filling ratio of ≥90% was set as a filter to single out the pixels in which rice paddies were almost completely covered.

Estimations of Phenological Metrics and Daily LAI Time Series in Paddy Fields
TIMESAT is a program that uses sophisticated statistical techniques to smooth satellite time-series data and retrieve key phenological parameters (http://www.nateko.lu.se/ TIMESAT/timesat.asp, access date: 15 September 2021) [34,35], which has been integrated in the processing of MODIS data into a phenology product (MOD09PHN and MOD15PHN) by the North American Carbon Program [36]. Boschetti et al. [37] compared the LGP and SoS derived from the TIMESAT and field census, respectively, by rice farmers, indicating a good correlation (R 2 = 0.92). Meanwhile, Son et al. [38] proved that rice sowing and harvesting dates could be exactly retrieved from MODIS VIs using the double logistic-based method. In line with these studies, the built-in noise reduction techniques in the TIMESAT program were used in this study to process time series of MODIS VIs, i.e., the spike method of median filtering, envelope iterations, minimum of 0.25, 0.1, and 0.05 in NDVI, EVI, and OSAVI, respectively, the fitting method based on the double logistic-based model, and the start of season method (i.e., STL trend). The study of forest ecosystems with remote sensing by Wu et al. [39] suggests that a better estimation of LGP and SoS was achieved by taking multiple VIs rather than one into account. Hence, NDVI, EVI and OSAVI were used to estimate seasonality parameters of paddy pixels. We defined the LGP per pixel, exceeding 160 days and less than 96 days, as outliers, because the modern semi-dwarf rice varieties in South Korea do not have LGP outside of the time range. Each seasonality index was determined using the average of three VIs estimates to lessen the impacts of random factors. In general, the changes in the LGP derived from the TIMESAT along with the growing season mean air temperature (gsT mean ) matched well with the regression line fitting the in-situ LGP-gsT mean data ( Figure S4a in Supplementary Materials SP2). The paired in-situ LGP and gsT mean datasets were collected from paddy rice fields [40]. A full agreement between the TIMESAT LGP and field observations suggests that the phenological metrics estimations of paddy rice using the TIMESAT were acceptable.
The daily LAI per-paddy pixel was generated according to the statistical relationships between VIs and LAI in paddy rice fields. Four-year measurements of ground surface LAI and VIs in three paddy sites in South Korea enabled us to perform a robust statistical correlation between the LAI and VI of interest. Three VIs (i.e., NDVI, EVI, and OSAVI), commonly used to retrieve LAI, were selected for comparison with field LAI observations. Results revealed the best fitting curve for the ln(LAI)-ln(OSAVI) correlation (R 2 = 0.72), followed by the ln(LAI)-ln(EVI) correlation (R 2 = 0.68), and the ln(LAI)-ln(NDVI) correlation (R 2 = 0.61) ( Figure S5 in Supplementary Materials SP3). Importantly, the variations of LAI > 1.0 m 2 m −2 (i.e., log(LAI) = 0) could be better explained by OSAVI. The same phenomenon for the raw data of LAI and OSAVI was also evident (The intent of Figure S5c). Therefore, seasonal courses of daily LAI for each paddy pixel in each year were estimated through the OSAVI time series. Beck et al. [41] demonstrated that the double logistic function can more accurately show seasonal changes of VIs by effectively removing outliers. Gonsamo et al. [42] and Wu et al. [39] examined a useful application of the double logistic model for noise reduction using remote sensing time series data. In our analysis, the double logistic model was fitted to time series of MODIS VIs data for each rice paddy and year, thereby yielding daily VIs as followings: where x 1 and x 3 are the locations of the left (the increasing phase) and right (the decreasing phase) inflection points respectively, and x 2 and x 4 determine the rates of changes at these points. Daily VIs values for each paddy field and year were then converted to LAI values based on the correlation between OSAVI and LAI.

Evaluation and Application of the Pixel-Based RS-PM Model for the ET Estimation in Paddy Fields
Daily ET in each paddy pixel was estimated by the pixel-based RS-PM model. The core principle for ET estimation in the pixel-based RS-PM model is the Penman-Monteithtype ET model modified by the FAO (i.e., FAO-56 PM model) [43]. It estimates ET based on a reference crop evapotranspiration (ET o ) multiplied by the sum of the transpiration coefficient (K cb ) and the evaporation coefficient (K e ) of the crop of interest. Our previous study proposed the pixel-based RS-PM model that integrated the FAO-56 PM model and the close-range unmanned aerial vehicle (UAV) based multispectral remote sensing for paddy rice grown at four fertilizer application rates in Gwangju, South Korea [12], thereby developing parameter retrieval formulas for K cb and K e . Detailed descriptions of the pixel-based RS-PM model structure with K cb and K e retrieval formulas can be found in the Supplementary Materials SP3.
To extract the spatial information on rice field ET over a large space with high simulation accuracy, the optimum configuration of key parameters for the FAO-56 PM model was necessary. As suggested by Allen et al. [43], the FAO-56 PM model could be performed well at the point and regional scales with the reasonable parameterization of independent weather and physiological variables. Sensitivity analyses of three important physiological parameters (i.e., aerodynamic resistance-r a , bulk surface resistance-r s , and K cb ) to estimate the daily ET estimate were performed in theory before field evaluations ( Table 1). The r a and r s describe the resistance of evaporating and transpiring surfaces to water vapor transport through the reference heights. The variation of r a from 20 s m −1 to 120 s m −2 caused a slight decline of daily ET by 10.02% (ranging from 4.09 mm m −2 day −1 to 3.68 mm m −2 day −1 ) under high-level radiation and the daily mean air temperature of 17.5 • C. The extents of ET constraint by escalating r a were 9.04%, 8.39% and 7.71% under varying daily air temperatures of 22.5 • C, 27.5 • C, and 32.5 • C, respectively. Results implied that configuring r a according to field conditions would not yield large errors for daily ET estimates. The varying r s values that affect ET through canopy transpiration caused a decline in daily ET at the same extent as increasing the r a . Significant impacts on daily ET estimates were evident in scenario 3 of Table 1, wherein varying K cb values, ranging from 0.58 to 0.30 caused a large decline in ET by almost 50%.
For the determination of the robust r s that is reliable for paddy rice ET estimation throughout South Korea, r s values increased from 20, 50, and 80 to 120 s m −1 , and the varying r s values fell within the ranges reported for paddy rice. Varying r a values were calculated on the basis of the seasonal crop height, according to the FAO-56 PM model instruction. The height of the rice crop was calculated using a relationship between LAI and plant height reported by Confalonieri et al. [44].  Table S2 in Supplementary Materials SP4). A time series of NDVI from 2002 to 2010 was obtained from gridded L3G (level3) composite data at 250 m spatial resolution (MODIS Terra Surface Reflectance products, MOD13Q1). The daily NDVI of each site was estimated through Equation (5) and used to evaluate seasonal K cb and K e . Comparisons between ET predictions and observations showed that the greater the r s , the smaller the root mean square error (RMSE). The FAO-56 PM model with a constant r s value of 120 s m −1 and the constant r a (called as FAO-56 PM rs = 120 ) had the best accuracy, with R 2 = 0.91 and RMSE = 0.50 mm m −2 day −1 in the Haean, R 2 = 0.71 and RMSE = 0.60 mm m −2 day −1 in the Mase, and R 2 = 0.64 and RMSE = 0.60 mm m −2 day −1 in the El Saler-Sueca. Relatively poor performances were observed when considering the variable r a and constant r s (120 s m −1 ) (termed FAO-56 PM 120,ra ). The RMSE of the FAO-56 PM rs = 120 model was 0.56 mm m −2 day −1 , accounting for around 17.5% of the mean daily ET (~3.2 mm m −2 day −1 ). As demonstrated in Table 1, we found that r s estimates would cause ET simulation errors of less than 10%. We found that estimation errors due to uncertainties in other model variables including r a would be lower than 10% as well; this is in agreement with sensitivity analyses of model parameters reported by our previous study [12]. Daily ET estimates in the Haean, Mase, and El Saler-Sueca over six years were compared favorably with EC observations (a in Supplementary Materials SP5). It was reasonable to infer that the pixel-based RS-PM model with optimization of parameter, including r s and r a enabled the efficient capture of the overall variation of rice field ET across South Korea. Table 1. The sensitivity analysis indicated the effects of parameterization themes of different models on daily ET estimates (mm m −2 day −1 ). The background environment in scenario 1 included incident solar radiation (R s ↓) of 19.0 MJ m −2 day −1 , daily relative humidity (RH) of 70% and the flooded bare land. Other model parameters were set at optimal levels except for the one used for the specific scenario analysis. In scenario 2, to evaluate the effects of canopy resistance (r s , s m −1 ), R s ↓ was set at 19.0 MJ m −2 day −1 , with daily RH of 70%, surface albedo of 0.18, and leaf area index (LAI) of 4.0 m 2 m −2 . In scenario 3, there were the same climate variables as in scenario 1 but with variying K cb values. Daily R s ↓ (MJ m −2 day −1 ) was estimated with data on weighted slopes, surface, and atmospheric moisture data by using multispectral imagery from the MI (Meteorological Imager) on COMS (Communication, Ocean and Meteorological Satellite), along with VEGETATION sensor data from SPOT (Système Pour l'Observation de la Terre) based on the Kawamura physical model [45]. The estimated daily R s ↓ was validated by observations made from 37 ground pyranometers, and the results showed a high correspondence between predictions and observations under all sky conditions [45,46]. Daily maximum and minimum air temperatures (T max , T min , respectively, • C) and the relative humidity (RH, %) obtained from the Korea Local Analysis and Prediction System (KLAPS) were used. The KLAPS, which constitutes the data collection and analysis modules, is designed to predict the weather conditions of the Korea Peninsula with a spatial resolution of 5.0 km at fine timescales. Weather products at high-resolution (1.5 km) from reanalysis data at 2.0 m height were generated by the KLAPS based on its analysis scheme using 470 automated weather systems throughout the Korea Peninsula [46,47]. Daily T s at 2.0 m height was the mean of daily T max and T min . Daily vapor pressure deficit (VPD, kPa) is a product of T s and RH. The same climate data as in our study were used for the estimation of ET and yield in paddy rice in South Korea by Yeom et al. [45] and Jeong et al. [46]. Wind speed grid data (m s −1 ) at a spatial resolution of 0.25 • and 1-day time intervals in South Korea were accessed from the Global Land Data Assimilation System Version 2 (GLDAS-2.0, Li et al., 2018) and downscaled with high-resolution covariates of advanced spaceborne thermal emission and reflection radiometer (ASTER) global digital elevation model (GDEM), slope, and aspect attributes through comprehensive statistical analyses using the thin plate smoothing splines and machine learning assembling (MACHISPLIN package in R software). The MACHISPLIN statistical technique was used to downscale the temperature and relative humidity products to match them spatially with MODIS VI products. Methods for daily R n estimation in this study were aligned with the FAO-56PM model. constant rs value of 120 s m -1 and the constant ra (called as FAO-56 PMrs = 120) had the best accuracy, with R 2 = 0.91 and RMSE = 0.50 mm m -2 day -1 in the Haean, R 2 = 0.71 and RMSE = 0.60 mm m -2 day -1 in the Mase, and R 2 = 0.64 and RMSE = 0.60 mm m -2 day -1 in the El Saler-Sueca. Relatively poor performances were observed when considering the variable ra and constant rs (120 s m -1 ) (termed FAO-56 PM120,ra). The RMSE of the FAO-56 PMrs=120 model was 0.56 mm m −2 day −1 , accounting for around 17.5% of the mean daily ET (~3.2 mm m -2 day -1 ). As demonstrated in Table 1, we found that rs estimates would cause ET simulation errors of less than 10%. We found that estimation errors due to uncertainties in other model variables including ra would be lower than 10% as well; this is in agreement with sensitivity analyses of model parameters reported by our previous study [12]. Daily ET estimates in the Haean, Mase, and El Saler-Sueca over six years were compared favorably with EC observations (Figures S6 and S7a in Supplementary Materials SP5). It was reasonable to infer that the pixel-based RS-PM model with optimization of parameter, including rs andra enabled the efficient capture of the overall variation of rice field ET across South Korea. Daily Rs↓ (MJ m -2 day -1 ) was estimated with data on weighted slopes, surface, and atmospheric moisture data by using multispectral imagery from the MI (Meteorological Imager) on COMS (Communication, Ocean and Meteorological Satellite), along with VEGETATION sensor data from SPOT (Système Pour l'Observation de la Terre) based on the Kawamura physical model [45]. The estimated daily Rs↓ was validated by observations made from 37 ground pyranometers, and the results showed a high correspondence between predictions and observations under allsky conditions [45,46]. Daily maximum and minimum air temperatures (Tmax, Tmin, respectively, °C) and the relative humidity (RH, %) obtained from the Korea Local Analysis and Prediction System (KLAPS) were used. The KLAPS, which constitutes the data collection and analysis modules, is designed to predict the weather conditions of the Korea Peninsula with a spatial resolution of 5.0 km at fine timescales. Weather products at high-resolution (1.5 km) from reanalysis data at 2.0 m height were generated by the KLAPS based on its analysis scheme using 470 automated weather systems throughout the Korea Peninsula [46,47]. Daily Ts at 2.0 m height was the mean of daily Tmax and Tmin. Daily vapor pressure deficit (VPD, kPa) is a product of Ts and RH. The same climate data as in our study were used for the estimation of ET

The Monthly Estimation of ∆T s between Paddy Field and Temperate Forest
α is an important component of surface energy balance, explaining the warming and cooling effects of paddy field. Daily values of α in paddy rice and temperate forest in South Korea were retrieved from MODIS MCD43A3 products. Daily T s and α values were aggregated into a value for the month (mean). In line with methods reported by Li et al. [3], the impacts of varying FFTD sat values of paddy field on monthly T s were expressed as the monthly T s difference of paddy rice minus that of the temperate forest.
Positive or negative ∆T s indicates a warming or cooling effect of paddy rice, respectively. Intensive paddy rice cultivation is primarily done in flat plains of the western part of South Korea (e.g., the Charyeong rice plain in Figure S1) as well as in mountainous areas. Therefore, we chose the Charyeong rice plain and surrounding forest ecosystems as two comparative groups (the rectangle extent: LL 35 • 8 N, 126 • 2 E, UR 37 • 3 N, 127 • 7 E;area: 1.8 × 10 6 ha −1 ; elevation range: 1 to 543 m; 337 forest pixels and 210 rice pixels in ). The average value of monthly ∆T s over four years is presented. Paddy rice is commonly raised on flat plains and forests at a wide range of elevations and slopes. Topography-induced differences in daily T s may enlarge or narrow ∆T s between paddy field and temperate forest. To correct the topography-induced biases in ∆T s estimation, we performed an elevation-aspect-slope adjustment by subtracting the elevation/aspect/slope-induced ∆T s from the original value, which was in line with the topography correction method proposed by Li et al. [3]. ∆α was the product of paddy field α minus that of the temperate forest. Spatial analysis indicated that cooling or warming effects on paddy rice were independent of the spatial extent selected for comparisons in South Korea.
Structural equation modeling (SEM) is a label for a diverse set of methods used in experimental and observational research, implying a structure for the covariances between the observed variables [48]. It was applied here for assessing relative importance of biophysical factors to regional variations of paddy field ET. All statistical analyses including SEM and coefficient of variation (CV) were analyzed by the R software (R Core Team, 2020) [49]. Spatial data processing was made by the ENVI-IDL 5.3 software (Exelis Visual Information Solutions Inc., Boulder, CO, USA).

Spatial Variations of ET in Paddy Fields, South Korea
The contribution of various biophysical factors to spatial variations of PF-ET and the rice-growing season total ET (RGS-ET, the sum of daily ET from the SoS to the EGS) is disentangled in this and the next paragraphs. The SEM was applied to evaluate the magnitude and importance of the causal connection between the RF-ET and nine biophysical factors (Figure 3), including RGS-ET, biomass production (estimation methods for biomass production referred to Supplementary Materials SP3), daily mean net radiation from the FFTD sat to the EGS, FFTD sat , daily mean and minimum air temperatures of the same period as the net radiation, SoS, LGP, the maximum LAI of the rice-growing season, and daily mean vapor pressure deficit during the same period. Explanatory power by the error item (PF-ET variations that could not be explained by the nine variables) was 0.024, suggesting that almost all spatial variations of PF-ET, reaching up to 98%, had been explained by these nine variables. The determination coefficient (DC), a product of the squared path coefficient, is a powerful indicator of the relative strength of an independent variable in determining the target variable. The DC value of ≥0.1 conventionally indicates a stronger correlation [48]. The DCs for the RGS-ET, SoS, and FFTD sat were 0.61, 0.12, and 0.18, respectively, remarkably greater than other values. A negative path coefficient between the PF-ET and FFTD sat was found, indicating that advancing the FFTD sat resulted in a greater PF-ET value. Positive signs for PF-ET and RGS-ET and for PF-ET and SoS suggested that the reductions in values of the two variables led to lower PF-ET value. The RGS-ET, FFTD sat , and SoS variables were critical (listed in order of importance) for measuring spatial variations of PF-ET, which was in agreement with the H1.
Spatial variations of RGS-ET could be well explained by eight biophysical factors (i.e., biomass production, daily mean net radiation during the rice-growing seasons-gsR n , daily mean minimum air temperature during the rice growing seasons-gsT min , SoS, daily mean air temperature during the rice growing seasons-gsT mean , LGP, growing season LAI max -gsLAI max , and growing season mean VPD-gsVPD), with an explanatory power up to 0.97. The DCs for the LGP, gsVPD, gsT mean , gsT min , and gsLAI max were 0.25, 0.13, 0.10, 0.11 and 0.10, respectively (Figure 3), suggesting that the LGP, gsVPD, gsT mean , gsT min and gsLAI max were the critical biophysical factors responsible for regional variations of RGS-ET.
The spatial CVs of PF-ET and RGS-ET and main biophysical factors causing differences in the CV between PF-ET and RGS-ET were quantified and explained in this paragraph. Figure 4 shows the regional variations of PF-ET and RGS-ET from 2011 to 2014. Clearly, the annual CVs of PF-ET tended to be conservative and constant at 0.09 throughout the four years. The annual CV of RGS-ET fluctuated from 0.11 to 0.19. Greater spatial variations of RGS-ET by 54.28% than those of PF-ET were observed. PF-ET was the summation of the rice-growing season total ET (i.e., RGS-ET) and the non-growing season ET (i.e., the period from the FFTD sat to the SoS). Results suggest that the field-to-field differences in the FFTD sat did not contribute to spatial variations of PF-ET. Nevertheless, the observed site-to-site differences in FFTD sat indeed lessen the expanding effects of RGS-ET spatial variations on PF-ET, as the spatial variations of RGS-ET were amplified by climatic and phenological factors, while the main regulating factors for PF-ET were FFTD sat and SoS.
by these nine variables. The determination coefficient (DC), a product of the squared path coefficient, is a powerful indicator of the relative strength of an independent variable in determining the target variable. The DC value of ≥ 0.1 conventionally indicates a stronger correlation [48]. The DCs for the RGS-ET, SoS, and FFTDsat were 0.61, 0.12, and 0.18, respectively, remarkably greater than other values. A negative path coefficient between the PF-ETand FFTDsat was found, indicating that advancing the FFTDsat resulted in a greater PF-ETvalue. Positive signs for PF-ETand RGS-ET and for PF-ETand SoS suggested that the reductions in values of the two variables led to lower PF-ETvalue. The RGS-ET, FFTDsat, and SoS variables were critical (listed in order of importance)for measuring spatial variations of PF-ET, which was in agreement with the H1.  (FFTDsat, DOY), the length of the growing period (LGP, days), the mean vapor pressure deficit (gsVPD, kPa) during the growing season, the daily mean air temperature (gsTmean, °C) during the growing season, the mean minimum air temperature (gsTmin, °C) during the growing season, the maximum leaf area index (gsLAImax, m 2 m -2 ) during the growing season, the daily net radiation over the rice-growing season (dRn, MJ m -2 day -1 ), the daily VPD(dVPD, kPa), the daily Tmin (dTmin, °C), the daily LAI (dLAI, m 2 m -2 ), the daily relative humidity (dRH, %) and the daily precipitation (dP, mm). Numerical value attached to the arrow symbol was path coefficient. Green (blue) sign indicated a positive (negative) correlation. DOY: the day of year. Spatial variations of RGS-ET could be well explained by eight biophysical factors (i.e., biomass production, daily mean net radiation during the rice-growing seasons-gsRn, daily mean minimum air temperature during the rice growing seasons-gsTmin, SoS, daily mean air temperature during the rice growing seasons-gsTmean, LGP, growing season LAImax-gsLAImax, and growing season mean VPD-gsVPD), with an explanatory power up to 0.97. The DCs for the LGP, gsVPD, gsTmean, gsTmin, and gsLAImax were 0.25, 0.13, 0.10, 0.11 and 0.10, respectively (Figure 3), suggesting that the LGP, gsVPD, gsTmean, gsTmin and gsLAImax were the critical biophysical factors responsible for regional variations of RGS-ET.
The spatial CVs of PF-ETand RGS-ETand main biophysical factors causing differences in the CV between PF-ETand RGS-ETwere quantified and explained in this paragraph. Figure 4 shows the regional variations of PF-ETand RGS-ETfrom 2011 to 2014. Clearly, the annual CVs of PF-ETtended to be conservative and constant at 0.09 throughout the four years. The annual CV of RGS-ETfluctuated from 0.11 to 0.19. Greater spatial Figure 3. Apath diagram presenting significant relationships and the causal direction for the paddy field total ET (PF-ET, the sum of daily ET from the FFTD sat to the EGS, mm m −2 ), and the rice-growing season total ET (RGS-ET, the sum of daily ET from the SoS to the EGS, mm m −2 ), with various biophysical factors, including the growing season onset of rapid canopy expansion (SoS, DOY), the first date of field flooding and transplanting identified by satellite observations (FFTD sat , DOY), the length of the growing period (LGP, days), the mean vapor pressure deficit (gsVPD, kPa) during the growing season, the daily mean air temperature (gsT mean , • C) during the growing season, the mean minimum air temperature (gsT min , • C) during the growing season, the maximum leaf area index (gsLAI max , m 2 m −2 ) during the growing season, the daily net radiation over the rice-growing season (dR n , MJ m −2 day −1 ), the daily VPD(dVPD, kPa), the daily T min (dT min , • C), the daily LAI (dLAI, m 2 m −2 ), the daily relative humidity (dRH, %) and the daily precipitation (dP, mm). Numerical value attached to the arrow symbol was path coefficient. Green (blue) sign indicated a positive (negative) correlation. DOY: the day of year.
No significant effects of spatial changes in FFTD sat on the LGP were observed ( Figure S4b), implying that the initial flooding activities would not bring about positive effects on RGS-ET through their impacts on LGP.

Seasonal Variations of ET in Paddy Rice Fields, South Korea
The contribution of various biophysical factors to seasonal variations of ET in regional paddy rice fields was disentangled in this paragraph. The analysis of the SEM model for daily ET during the rice-growing seasons indicated that seasonal variations were primarily related to dLAI, dVPD, dT min and dR n (Figure 3), with the relative importance scores in the following order: dR n > dVPD > dT min > dLAI. To take account of the potential effects of other meteorological factors such as soil temperature and daily rainfall that were not available for the gridded ET products, ground-surface measurements in the Haean, Mase, and El Saler-Sueca were supplemented. New SEM analysis based on in-situ observations showed that the daily ET was mainly controlled by dVPD and dR n , followed by dT min and dLAI. The latter two parameters made smaller contributions to daily ET from the FFTD to the EGS. Results derived from the RS-PM model were comparable to those of SEM analysis based on EC observations. Additionally, the significantly negative and positive signs for the dVPD-dRH and dRH-dT min relationships, respectively, were found. The causal coefficient for the dT min -dVPD relationship was much smaller (data are not presented). Hence, the positive effects of dT min on daily ET resulted by its influences on dRH thereby dVPD. We speculate that the significant effects of air temperature on spatial variations of RGS-ET were produced through their controls on dVPD. mation of the rice-growing season total ET (i.e., RGS-ET) and the non-growing season ET (i.e., the period from the FFTDsat to the SoS). Results suggest that the field-to-field differences in the FFTDsat did not contribute to spatial variations of PF-ET. Nevertheless, the observed site-to-site differences in FFTDsat indeed lessen the expanding effects of RGS-ETspatial variations on PF-ET,as the spatial variations of RGS-ETwere amplified by climatic and phenological factors, while the main regulating factors for PF-ETwere FFTDsat and SoS.

Seasonal Variations of ET in Paddy Rice Fields, South Korea
The contribution of various biophysical factors to seasonal variations of ET in regional paddy rice fields was disentangled in this paragraph. The analysis of the SEM model for daily ET during the rice-growing seasons indicated that seasonal variations were primarily related to dLAI, dVPD, dTmin and dRn (Figure 3), with the relative importance scores in the following order: dRn>dVPD>dTmin>dLAI. To take account of the potential effects of other meteorological In this and the following paragraphs, daily ET time series from the FFTD sat to the EGS in paddy field with two kinds of FFTD sat and biophysical mechanisms causing differences in ET variation between the rice-growing and the non-growing seasons were quantified. The seasonal pattern of daily ET was generated by taking the average of seasonal ET value in all rice paddies across South Korea. Temporal changes of daily ET in 145 DOYfirst flooding paddies (145 FFTD sat ) were compared with those in 121 DOY-first flooding paddies (121 FFTD sat ) ( Figure 5). The maximum daily ET of~6.0 mm m −2 day −1 was recorded in the mid-season. Daily ET in two types of flooding fields fluctuated greatly from the FFTD sat to the EGS and throughout the four years. LAI rapidly increases after rice transplanting, arrived at the maximum values around mid-July (200 DOY), and decline thererafter. Seasonal daily ET change in 2012 and 2013 did not matched well with LAI development. Decoupling between daily ET and LAI in Figure 5d was consistent with EC observations in the Haean and Mase sites ( Figure S6a-d).
(145 FFTDsat) were compared with those in 121 DOY-first flooding paddies (121 FFTDsat) (Figure 5). The maximum daily ET of ~6.0 mm m -2 day −1 was recorded in the mid-season. Daily ET in two types of flooding fields fluctuated greatly from the FFTDsat to the EGS and throughout the four years. LAI rapidly increases after rice transplanting, arrived at the maximum values around mid-July (200 DOY), and decline thererafter. Seasonal daily ET change in 2012 and 2013 did not matched well with LAIdevelopment. Decoupling between daily ET and LAIin Figure 5d was consistent with EC observations in the Haean and Mase sites (Figures S6a-d).  SoS, as an indication of the rapid development of canopy leaf area, was intensively distributed between 160 DOY and 180 DOY ( Figure S8 in Supplementary Materials SP6). Hence, we chose SoS at the leftmost boundary of 160 DOY as a criterion to separate the ricegrowing season from the non-growing season. The non-growing season from the FFTD sat to 160 DOY was characterized by low LAI value, and thereby a large proportion of land surface was directly exposed to air. For seasonal variations of daily ET in the 121 FFTD sat paddies, the CV of ET during the non-growing season was 0.28 that was substantially lower (by 28.2%) than that during the rice-growing season ET (CV = 0.39) (the rice-growing season was from 160 DOY to the EGS) in 2011. The CVs of ET during the non-growing season in 2012, 2013, and 2014 were 0.23, 0.29, and 0.21, respectively, significantly lower than the corresponding CV of ET during the rice-growing season. Greater seasonal variations of ET during the rice-growing season compared to that during the non-growing season ET were repeatedly observed for the 145 FFTD sat paddies in the same years. The extending effects of greater spatial and seasonal variations of ET during the rice-growing season on the PF-ET were significantly minimized by advancing the FFTD sat .

Spatial Distribution Attributes of Paddy Field FFTD sat in South Korea
Spatial distribution of the FFTD sat over paddy rice fields of South Korea was mapped in Figure 6. Each 8-day increment was marked with a different color, with color gradients marking from the beginning of May (113 DOY) to the beginning of June (161 DOY). In the case of the field irrigation management in 2011, red, blue, and cyan dots were dominant color types and the other three were interspersed with wide space. Intensive rice paddy systems were found on the midwestern plains where the elevation was lower than 100 m. Some rice paddies had the FFTD sat scheduled on 129 DOY on midwestern plains in 2011, while the FFTD sat for most rice paddies in 2012 shifted to 113 DOY and 137 DOY, as the surface areas of red and cyan dots were greater in the same area. Apparently, red, blue, and cyan dots co-existed on the mid-west plains in 2013 and 2014, indicating the range of FFTD sat from late April to the end of May; this was in agreement with our previous research reports for the same area and year (please refer to Jeong et al., 2018) [46]. Meanwhile, FFTD sat frequency distribution in Figure 6 indicated that 90% of sampling paddy fields had the FFTD sat values ranged from 110 DOY to 140 DOY, suggesting that field flooding and transplanting campaigns were intensively initiated from the end of April to the end of May in most paddy fields in South Korea.   (Figure 7a). The 145 FFTD sat paddies had a higher α value (28.8%) than 113 FFTD sat fields in May, implying that advancing the FFTD sat could result in a low surface reflectance and more solar energy absorption. In June, the α value of the temperate forest increased to 0.143, while it remained at a lower level (37.1%) in the 113 FFTD sat fields, with a negative ∆α (Figure 7a). The α value of 145 FFTD sat fields experienced a dramatic decline (16.4%) from May to June, as the 145 FFTD sat fields were barren/fallow in May and changed to flooding paddies thereafter in June. Large differences in α value between 113 FFTD sat and 145 FFTD sat fields in May narrowed in June without a statistical significance at 0.05 level (p-value = 0.10). From July to September, α value reached its peak (0.168), while it did not show significant differences between 113 FFTD sat and 145 FFTD sat fields thereafter (p-value > 0.05). Intra-annual changes in the maximum/minimum air temperature difference (ΔTmax, ΔTmin, respectively) between paddy rice and the temperate forest are presented in Figure 7b,c. The monthly ΔTmax had a bottom-down bowl curve or V-shaped pattern with positive values at the start and end of the year and negative values from April to July. The lowest negative value (−1.25 °C) was recorded in May when flooding had started or finished in most paddy rice fields. The time node when the negative sign of ΔTmax shifted to positive was August. Although positive signs were found in August and September, the ΔTmax values in the two months were relatively lower than 0.25 °C. The ΔTmin value remained positive from January to the end of the year, with the highest values recorded in April and October (Figure 7c). Except for the winter, ΔTmin values were greater than 0.5 °C in the other three seasons. Taking the mean of ΔTmax and ΔTmin into consideration, seasonal changes in profiles of ΔTs were observed (Figure 7d), showing a shallow bottom-down bowl shape similar to that for the ΔTmax. However, negative ΔTs values were found in May and June, and a positive value slightly lower than 0.25 °C was observed in July, wherein ΔTmax was negative. It was obvious that the nighttime warming amplitudes in April and July outpaced the daytime cooling, resulting in the positive ΔTs values in two months.
This paragraph compared the ET values in paddy rice and the temperate forest. We cut outthe MOD16A2grid ET data for the 3 × 3 pixelssurrounding the ECtower of each site. Seasonal changing profiles of 8-day ET (the sum of daily ET for every eight days) derived from three data sources (ECtower, MOD16A2 aET-actual ET, and RS-PMestimations) are shown in Figure S10 in Supplementary Materials SP7. A close correspondence in 8-day ET between EC observations and RS-PMsimulations was observed. However, MOD16A2 aET products could not accurately regenerate 8-day ET of paddy rice in three sites, especially during the land preparation process and early growth stages of rice with distinct underestimations by two to three times (the shaded area in Figure S10). The daily Intra-annual changes in the maximum/minimum air temperature difference (∆T max , ∆T min , respectively) between paddy rice and the temperate forest are presented in Figure 7b,c. The monthly ∆T max had a bottom-down bowl curve or V-shaped pattern with positive values at the start and end of the year and negative values from April to July. The lowest negative value (−1.25 • C) was recorded in May when flooding had started or finished in most paddy rice fields. The time node when the negative sign of ∆T max shifted to positive was August. Although positive signs were found in August and September, the ∆T max values in the two months were relatively lower than 0.25 • C. The ∆T min value remained positive from January to the end of the year, with the highest values recorded in April and October (Figure 7c). Except for the winter, ∆T min values were greater than 0.5 • C in the other three seasons. Taking the mean of ∆T max and ∆T min into consideration, seasonal changes in profiles of ∆T s were observed (Figure 7d), showing a shallow bottom-down bowl shape similar to that for the ∆T max . However, negative ∆T s values were found in May and June, and a positive value slightly lower than 0.25 • C was observed in July, wherein ∆T max was negative. It was obvious that the nighttime warming amplitudes in April and July outpaced the daytime cooling, resulting in the positive ∆T s values in two months.
This paragraph compared the ET values in paddy rice and the temperate forest. We cut outthe MOD16A2grid ET data for the 3 × 3 pixels surrounding the EC tower of each site. Seasonal changing profiles of 8-day ET (the sum of daily ET for every eight days) derived from three data sources (EC tower, MOD16A2 aET-actual ET, and RS-PM estimations) are shown in Figure S10 in Supplementary Materials SP7. A close correspondence in 8-day ET between EC observations and RS-PM simulations was observed. However, MOD16A2 aET products could not accurately regenerate 8-day ET of paddy rice in three sites, especially during the land preparation process and early growth stages of rice with distinct underestimations by two to three times (the shaded area in Figure S10). The daily mean ET values of the forest land derived from the MOD16A2 products were 1.85 mm m −2 day −1 in May and 3.10 mm m −2 day −1 in rice paddies according to the RS-PM model and EC observations. The daily mean ET value of the temperate forest in August and September was 2.1 mm m −2 day −1 according to the MOD16A2 aET, which was higher than that of senescing rice plants, aligning with actual recognitions. Although MOD16A2 products failed to capture daily ET values of rice paddies, we confidently believe that the daily ET value of the temperate forest in May might be lower than that of flooded paddy fields, as the abundant water resources were available for rice paddies which facilitate the water loss through ET, and the forest was going through the greening stage during the spring with small green leaf surface area in canopies featured by low canopy transpiration rates. The partial flooding of paddy fields in April and the intensive flooding in May incur water loss at daytime through higher ET values, dissipating excess energy absorbed during the daytime. Relative higher ET fluxes of rice paddies in May and June caused the decline of daytime air temperature that exceeded night-time warming amplitudes, leading to strong cooling effects in May and June. The results of the present study support and interpret the H2 stating that the monthly T s value of paddy fields was not always higher than that of the temperate forest. Increasing ET capability in August, September and October by forest ecosystems narrowed the differences in ∆T max , which further increased the positive ∆T s value in the three months in the context of positive ∆T min .

Field-to-Field Variations of FFTD sat Among Paddy Fields
Paddy rice plains are primarily distributed in the western part of the Korean Peninsula with flat topography and a land slope of lower than 10%. Some small-scale paddy fields are located in the middle and eastern parts, where the forest is the dominant land-use type. Paddy rice areas per county quantified by the flooding single classification method had a strong correlation with national census data with R 2 of up to 0.85 ( Figure S11 in Supplementary Materials SP8). Differences in monthly surface albedo between 113 FFTD sat , 145 FFTD sat paddy fields, and the forest due to the initial flooding serve as important evidence to support the capability of the flooding single classification method for accurate identification of FFTD sat in paddy rice fields. Small paddy rice landholdings and their FFTD values were characterized by the flooding single classification method in paddy rice fields.
Small paddy rice landholdings in East and Southeast Asia commonly receive various irrigation and fertilizer management practices. The water requirement for land preparation is theoretically 150-200 mm, but can be as high as 650-900 mm when its duration is prolonged to 30 days [50]. Jeong et al. [18] mapped the spatial distribution of FFTD sat in paddy rice in South Korea but did not analyze the reasons. The land preparation process of paddy fields represented by the FFTD sat data started from the end of April until the beginning of June, during which the weather was moderately dry with monthly average precipitation of less than 80 mm [51], that makes it impossible to provide enough water for land preparation; thus, water sources except for precipitation are indispensable for rice land preparation. The spatial analysis for rice paddies and river/stream systems showed that approximately half of the rice paddies (52%) were located within a 3.0 km range centered on permanent river/stream systems ( Figure S12a in Supplementary Materials SP9); however, another half is distributed far alway from the permanent river/stream systems. The statistical analysis for FFTD sat and the linear distance from river systems indicated no significant correlation between them ( Figure S12b). Pearson's correlations for FFTD sat and the daily air temperature averaged over a window size of 15 days ahead of the FFTD sat indicated the positive correlations between the two with correlation coefficients greater than 0.7 for all time windows in 2011 and greater than 0.93 in 2013 (Table S3). Nevertheless, significant correlations for the FFTD sat and the daily air temperature of 15-day time windows were not repeatedly presented in 2014. There was no significant correlation between FFTD sat and the air temperature before the first flooding date. It may suggest that the spatial distance of rice fields from river/stream systems may not be the main factor affecting the time of rice land preparation.
The previous study reported the effect of the observed site-to-site variations on the transplanting date due to prior adverse weather conditions such as frost risk [19]. In northern temperate areas, the rice land preparation usually starts when the daily air temperature is stably higher than 12 • C [52]. This indicates that the first flooding of paddy rice can be initiated at any day of daily air temperature stably >12 • C. As suggested by Jeong et al. [18] and Xiao et al. [29], personal preference, the empirical knowledge of individual farmers, was the key factor influencing the FFTD. However, our results showed that FFTD did not have a strong and consistent relationship with the daily air temperature during the land preparation from 2011 to 2014 (Table S3). It suggested that the air temperature in earlier days may not be the key factor considered by farmers when taking the initial flooding into account.
According to the annual statistical report for 2018 by OECD, the general service support transfers, which were linked to measures setting out enabling conditions for primary agricultural sectors through the development of water infrastructure and private or public services, amounted to 3639 million USD; thus, ranking third among the thirty-five OECD member countries [53]. Approximately, more than 80% of Korea's paddy production is supported by over 18,000 agricultural reservoirs and 70,000 other irrigation facilities [20]. The well-maintained and accessible water infrastructure in the agricultural sector could protect cropland against water-deficit stress caused by prolonged dry periods, which might be one reason for the weak correlations between FFTD, the linear distance of rice fields from river/stream systems, and air temperature. Another reason might be the difference in the ability of seedling to tolerate cold [54]. They found the pattern of latitude-dependent variation in cold tolerance in cultivated rice. Cold-tolerant cultivars might be introduced in early spring, which could result in weak correlations between FFTD and the air temperature during the early days. Further investigation to explain why the first flooding time over space spanned a long period (i.e., from the beginning of May to the beginning of June) was useful.

Spatiotemporal Variations of ET in Paddy Fields
Several studies have been conducted on daily ET during the rice-growing season [5][6][7]9,[11][12][13][14][15]. Researchers have arrived at the consensus regarding the daily ET at the initial stage of rice growth with the low leaf area index, mainly depending on climatic factors. During the mid and late-season stages in a closed canopy with a high leaf area index, transpiration through the rice leaf and stem is the dominant component of daily ET. The daily variability of ET during the two stages is considered to be dependent on climatic factors, canopy structure, and stomatal conductance. Among meteorological factors, incident solar radiation was considered to be critical, and leaf area index had a significant regulating effect on daily ET when it kept increasing. Our results showed that the daily ET during the rice-growing season was strongly correlated with the net radiation and water vapor pressure deficit; these findings were in agreement with those of other studies. Evaporation, as the main component of ET during the early growing stage and the non-growing period (from the FFTD to the SoS), was comparable to and even higher than the ET at mid-growth stage when LAI reached its maximum value. Hence, there was no considerable difference in daily ET between the non-growing period and the mid-growth stage, resulting in a weak relationship between daily ET and LAI. Seasonal effects of LAI on daily ET during the rice-growing season were concealed by meteorological factors due to the prolonged non-growing period as a result of advancing the FFTD.
The spatial variation of RGS-ET was larger than that of the PF-ET by 54.28%. Siteto-site changes in RGS-ET were primarily interpreted by the differences in LGP, gsVPD, gsT mean , gsT min , and gsLAI max .
LGP had the greatest contribution to spatial changes of RGS-ET, followed by VPD and air temperature. This finding was consistent with the results of a recent study on the spatial distribution of ET during the mid-season by Han et al. [55], who reported the temperature as an important meteorological factor. Nevertheless, the attribution of these factors to spatial variations of PF-ET followed the order of importance as follows: FFTD > SoS. Delaying the SoS and advancing the FFTD increased the PF-ET. The water loss through the evaporation from the FFTD to the SoS (i.e., the non-growing season) was the dominant component of ET since rice canopy was relative short, and a large proportion of land surface was directly exposed to the air during the non-growing season. Advancing the FFTD resulted in greater PF-ET via a prolonged day length during the nongrowing season, which subsequently caused a decline in the mass proportion of RGS-ET to PF-ET. Temporal fluctuations of daily ET during the non-growing season were much smaller than those during the rice-growing season, as the relatively stable incident solar radiation was the main factor controlling daily ET during the non-growing season. Greater spatial changes in RGS-ET that were not observed for PF-ET, were alleviated by prolonging the non-growing season via advancing the FFTD. It was reasonable to infer that human activities (i.e., the varying FFTD) alleviated spatiotemporal changes in evapotranspiration, which were amplified by the climatic and phenological changes in paddy rice ecosystems.
Wang et al. [19] reported that the earlier shift of the transplanting date (−2.0 ± 4.8 days/decade) alone prolonged the LGP of early rice by 1.3 ± 5.5 days/decade. They defined LGP as the number of days from the transplanting date to the end of the growing season, different from that reported in ours and others' studies on the flux ecology [39] and global change ecology [56]. We found a negative correlation between gsT mean and LGP ( Figure S4a), namely for each one-degree increase in the daily air temperature, and a decrease in LGP by 3.4 days was observed. This was supported by the phenology growing degree days theory that the accumulation of effective heat units accelerates the phenological development and lessens the LGP [57]. The LGP in our research was not correlated with FFTD, partially due to the marginal effects of FFTD on spatial variations of the air temperature during the rice-growing season.
Bringing the FFTD towards early spring induced a great water loss via ET during the non-growing season. The earliest flooding date reported in our analysis was around 113 DOY, almost 40 days earlier than that of paddy fields, with the first flooding dates scheduled at the end of May (e.g., 150 DOY). Nevertheless, grain production was not significantly influenced by the site-to-site differences in FFTD sat ( Figure S13 in Supplementary Mat SP10). On average, the amount of water loss by ET was~3.7 mm m −2 day −1 , which induced an annual unproductive water loss of 9.53 × 10 6 tons in South Korea. We recommended that flooding the rice fields followed by seedling transplantation in mid or late May might be the controllability measures to minimize the unproductive water loss in South Korea.

Human-Induced Seasonal Perturbations of Cooling and Warming Effects in Paddy Fields
The paddy rice cropping system differs substantially from any other upland nonirrigated or irrigated croplands. The distribution characteristics of small rice paddies mosaic within agroforestry basins and the large-scale intensive cropping are believed to bring about substantial influences on seasonal and annual dynamics of surface energy budgets, thereby affecting surface climate. Non-radiative warming of annual surface temperature by 0.5 • C was observed as a result of the global assessment of the forest conversion to cropland (including paddy rice) in the mid-latitudes [2,3]. As compared to cropland, the forest land in low and mid-latitude areas was considered to reduce the land surface temperature more significantly than the crops. In our research, the effects of the flooding campaign in rice paddies on intra-annual changes of the monthly T s were determined.
We found the difference in the monthly mean air temperatures between paddy rice and the temperate forest of −0.02 • C in May, −0.07 • C in June, and 0.15 • C in July, which were comparable to those in the temperate forest, implying a significant cooling effect of the paddy fields within the three months. Weather station observations between the wetland and alfalfa land in California reported seasonal oscillations of ∆T s between 1.0 and −1.0 • C [58]. Using station observations, Christy et al. [59] found no significant changes in time series of air temperature during the past century in the California Central Valley, on account of a counteracting effect by an approximate 3 • C decrease in the maximum near-surface temperature and a 3 • C increase in the minimum temperature. The maximum daytime air temperature in rice paddies was reduced by 1.25 • C in May, while the minimum night-time air temperature increased by 1.23 • C. The remarkable daytime cooling and nighttime warming effects of rice paddies in May were also reported by Liu et al. [25] and Yu and Liu [26]. Significant daytime cooling effects were canceled out in summer and autumn from August to October. Daytime warming effects in paddy fields vs. the temperate forest in September were observed, likely due to the lower ET in senescing rice canopy and enhanced ET in the temperate forest. As suggested by Li et al. [3], the ongoing evapotranspiration at night-time in the temperate forest dissipated more energy that was absorbed during the daytime, causing an annual net cooling effect as compared to the cropland. We observed the year-round night-time warming in paddy fields compared to the temperate forest, which might be due to the ongoing evapotranspiration at night in the temperate forest.
The differences in land surface temperature are a result of changes in attributes of surface energy balance, including soil heat flux, surface albedo, roughness, and the energy partitioning between latent and sensible heat flux [22]. The soil heat flux was supposed to be neglected when it lasted over 24 h or for a longer time [43], especially in paddy fields [9,11]. Fallow lands showed substantial declines in the surface albedo and improvement of daily ET after field flooding in May and June. A reduction of 22% in surface albedo was observed in early flooded paddy fields. With lower albedo in May and June, paddy fields absorbed more shortwave radiation during the daytime, potentially leading to a warming effect. However, this net energy gain was offset by a greater latent heat loss in the form of higher ET in flooded paddy fields, resulting in a daytime cooling effect of up to 2.1 • C as compared with the non-flooded barren fields. Paddy fields with the standing water layer and an open water surface serve as a large-size water body with a high heat capacity. The high heat capacity allows the puddled fields to lose heat more slowly while still staying warm at night. If the available solar energy is more, the paddy fields receive more energy via the lower surface reflectance during the daytime, resulting in a stronger night-time warming effect. Finally, we concluded that the early flooding schemes in rice fields have gained cooling benefits in the late spring and early summer, primarily due to a significant decline in the daytime air temperature as a result of the enhanced latent heat fluxes through ET.

Conclusions
The spatial analysis technique based on earth observations and numerical analyses is a relatively new approach to explore the spatiotemporal patterns of vegetation ET at large spatial scales. In our research, a data fusion technique that integrated earth observations (i.e., the MODIS surface reflectance, albedo and daily ground surface climate datasets) and numerical algorithms (i.e., the flooding single paddy rice classification method and the RS-PM model) provided us with new insights. A large number of studies have found seasonal cooling effects of paddy rice in Northeast and Southeast China, especially in May and June [25][26][27][28]. They reported that changes in cropland structure from the rain-fed farmland to the paddy field could profoundly modify land-surface thermal processes. Significant cooling effects of paddy rice in May and June were also found in South Korea in our research. Results in our study quantified the relative importance and contributions of biophysical factors to regional variations of ET and seasonal effects of temperature fluctuations in paddy fields, which advances our understandings of the impacts of land use and management change on regional climate change.
Results revealed that the spatial variations of PF-ET were substantially lower than those of RGS-ET. SoS and FFTD sat were the most important factors responsible for the spatial variations of PF-ET, while LGP, gsVPD, gsT mean , gsT min , and gsLAI max were predominant factors controlling spatial variations of RGS-ET. Spatial variations of PF-ET that were not as high as those of RGS-ET were mainly attributed to bringing the FFTD sat forward, as it prolonged the non-growing season during which the variation of daily ET was much lower than that during the rice-growing season. There was a significant cooling effect by early flooding in paddy fields compared to the temperate forest, primarily due to the large decline in maximum daytime air temperature. The magnitude of the monthly cooling effect in paddy fields was equal to that of the temperature forest, especially from May to June. Human activity-induced cooling effects in paddy fields showed equal importance as compared to the temperate forest, which could largely mask the ongoing global warming in spring. Further field studies designed to evaluate the social and biophysical factors influencing the site-to-site differences in FFTD in paddy fields across space would be useful, especially for deepening our understandings of spatiotemporal variations of ET in regional paddy fields and their impacts on climate change.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/rs13193992/s1, Figure S1: 5-m resolution raster image of main land-cover types in 2011, South Korea; Figure S2: Annual changes of paddy rice land use in main rice plain of Western Korea; Figure S3: Comparisons between map of MODIS potential paddy pixels and map of paddy rice land use and cover at 5-m spatial resolution; Figure S4: (a) Linear correlation between the length of growing period (LGP, days) and mean air temperature of the rice-growing season (gsT mean , • C) and (b) effects of the field first flooding and transplanting date detected by satellite observations (FFTD sat ) on the LGP; Statistical correlations between (a) ln(NDVI) and ln(LAI), (b) ln(EVI) and ln(LAI), (c) ln(OSAVI) and ln(LAI) in paddy fields in South Korea; Table S1: Selected equations involved in the GRAMI-rice model and its description; Table S2: Climate conditions at Haean (South Korea), Mase (Japan), and El Saler-Sueca (Spain); Figure S6: Seasonal changing courses of daily ET derived from EC systems installed in Haean, Mase, El Saler-Suecasites and the optimized RS-PM rs=120 model simulation for each site and year; Figure S7: Comparisons between predictions and observations of daily ET and grain yield in Japan_Mase, Korea_Haean, Spain_El Saler-Sueca, and Korean_Gwangju; Figure S8: Histogram charts of growing season onset (SoS) in paddy fields. Each seasonal metric corresponding to one flood date was the mean over 2011-2014; Figure S9: Surface albedo (α) responses to land cover and management change; Figure S10: Comparisons of 8-day ET of paddy fields in Haean, Mase, and El Saler-Sueca between EC observations, MOD16A2 products (aET and pET: actual ET and potential ET), and RS-PM rs=120 simulations; Figure S11: Comparison in paddy rice area per county between national statistic census and estimations by the flooding single-based classification method; Table S3: Pearson correlations for the first flooding and transplanting date estimated by satellite observations (FFTD sat ) and daily air temperature over different time windows; Figure S12: (a) Percentage of rice paddy with buffer distance away from permanent river/stream systems and (b) spatial correlation for the first flooding and transplanting date (FFTD sat ) and linear distance away from river/stream systems; Figure S13: Effects of the first flooding and transplanting date (FFTD sat ) on grain yield, PF-ET and WUE agro in 2011, 2012, 2013 and 2014, South Korea.