Estimating the Gross Primary Production and Evapotranspiration of Rice Paddy Fields in the Sub-Tropical Region of China Using a Remotely-Sensed Based Water-Carbon Coupled Model

: Rice serves as the staple food for over 50% of the global population. Remotely-sensed based estimation of the gross primary production (GPP) and evapotranspiration (ET) of rice paddy ﬁelds is essential to assess global food security. In this study, we tested the application of a recently proposed remotely-sensed based water-carbon coupled model (PML-V2) in the lower reaches of the Poyang Lake plain, which is one of the nine production bases for crops in China. Evaluation using the eddy covariance measurements showed that, after parameter localization, the model reproduced the seasonal variations of GPP and ET for both the early rice and the late rice. The model performed reasonably well in the validation period because the key parameters (e.g., the quantum efﬁciency and the stomatal conductance coefﬁcient) exhibited predictable seasonal variations. At the regional scale, the spatial distribution in multi-year GPP of rice (1365 ± 326 gCm − 2 year − 1 ) can be explained by the vegetation cover fraction (R 2 > 0.9); in comparison, the multi-year ET (1003 ± 65 mm/year) exhibits smaller spatial variations due to the high evaporation rate of the saturated soil surface of paddy ﬁelds. The water use efﬁciency of rice in this region varies around 1.35 gC/kgH 2 O with a standard deviation of 0.30. Our study shows that GPP and ET of rice can be estimated by remote sensing models without detailed crop management information, which is usually unavailable at regional scales.


Introduction
The fast-growing population is driving a fast-growing need in the food supply globally. In 2050, the increase of new cropland is projected to be about 1 billion ha [1], and the food production needs to be 70% higher to meet the caloric requirements of the global population [2]. Among all the staple foods, rice is a particularly important one because it feeds over 50% of the population worldwide [3]. Therefore, accurate estimation of rice yield and its water use, especially in the Asian monsoon region in which over 90% of the rice grows [4], is essential to guarantee global food security and sustainable development.
With the development of remote sensing techniques, numerous crop models [5] and statistical models [6] have assimilated the remotely sensed vegetation information to estimate regional crop yields [7]. For example, Sasai et al. [8] estimated the spatial and temporal variations in the rice productions in Japan with a diagnostic biosphere model and the Moderate Resolution Imaging Spectroradiometer (MODIS) leaf area index (LAI) product. The point-scale evaluation showed that, after introducing the harvest and soil oxidation-reduction processes, the root mean square errors were reduced to 27.31 gCm −2 month −1 and 36.33 gCm −2 month −1 for gross primary production (GPP) and net ecosystem production (NEP), respectively [8]. Other diagnostic models, such as the light use efficiency (LUE) model [9] and the light-response curve model [10] that were proposed to estimate ecosystem GPP without parameterizing the crop management information and the process at the vegetation-soil interface, were also widely used in the regional estimation of crop yields, for example, for the maize ecosystem [11], the wheat-maize rotating system [12], and rice paddy fields [13,14].
Irrigated croplands consume most of the precipitation by evapotranspiration (ET) [15] and therefore have strong feedback to the terrestrial hydrology and the near-surface atmosphere [16,17]. Compared to other crops and natural vegetation, the ET of the paddy field is generally much higher [18]. For example, Zhao et al. [19] reported that the ratio between ET and net radiation increases from 0.6 to 0.8 when a marshland was turned into a rice paddy field in the northeast region of China. Similarly, ET of rice takes up 63-67% of the total ET during the growing season for a rice-wheat rotation site in Nanjing, Eastern China [20]. Therefore, accurate estimation of regional ET for rice is essential for water management and sustainable agricultural development. The ET of rice can be estimated using the Priestley-Taylor (PT) method [21] or the Penman-Monteith model [22,23]. For example, Qiu et al. [24] examined the sensitivity of rice paddy ET to varying patterns of warming using a modified PT model in the East Asian monsoon regions. Xu et al. [25] estimated rice ET under water-saving irrigation by the PM model after calibrating the canopy resistance model parameters, in which LAI is used to scale up the stomatal conductance for water vapor from the leaf scale to the canopy scale.
To correctly assess the water use efficiency of crops, the rates of photosynthesis and transpiration of vegetation are usually estimated in a coupled way [26,27] because both carbon dioxide and water vapor go through leaves via the stomata. In addition, the stomatal conductance is sensitive to the photosynthesis rate [28]. Recently, Zhang et al. [29] and Gan et al. [30] proposed a routinely applicable water-carbon coupling model (PML-V2), which combines the carbon and vapor fluxes through the response of the canopy conductance to the photosynthesis rate. The PML-V2 model estimates GPP and ET using the light-response curve method and the PM model, respectively. Although the PML-V2 model has been tested at 95 EC sites globally [29], the usage of the model is seldomly evaluated for rice paddy fields. In this paper, based on eddy covariance (EC) data from a double-cropping rice paddy field in the lower reaches of the Poyang Lake basin, China, we intend to (1) test the usage of the PML-V2 model for a double-cropping rice paddy ecosystem in this area, (2) study the differences in key parameters of the model (light-responses, canopy conductances, and water-carbon coupling characteristics) between early rice and late rice, and (3) upscale the in situ GPP, ET, and water use efficiency (WUE) to the regional scale in the lower reaches in the Poyang Lake plain using the PML-V2 model.

Study Site
Our study area is located in the lower reaches of the Poyang Lake plain (Figure 1), which is one of the nine production bases for crops in China. Characterized as the humid subtropical climate, the annual mean air temperature is 17.5 • C, and the mean annual precipitation is 1635.9 mm for the Poyang Lake basin for 1960-2010 [31,32]. The major rain season of this area is April to June, and precipitation decreases sharply after July [32]. In this paper, we applied the model at a regional scale (28 • -29 • N, 115.7-116.7 • E) (Figure 1), where rice paddy fields occupy over 45% of the total land coverage, as shown in Figure 1. In addition, the selected area mainly contains the rice pixels that are near Poyang Lake, and the pixels in the Ganjiang and Fuhe watersheds, which are two major sub-watersheds of the Poyang Lake basin. Rice is the main crop in the Poyang Lake basin, where dual cropping systems are prevalent. The early rice is generally growing from late March to middle July, and the late rice is generally growing from middle July to late October. After the late rice is harvested, herbaceous plants, such as the Chinese milk vetch, are usually grown to maintain the fertility of the land. The typical cropping system and management measures [33] in this area are summarized in Table 1. In this paper, we will evaluate the performances of the PML-V2 model in simulating the GPP and ET at the Xiangtang station (28 • 26 25 N, 116 • 00 01 E), which belongs to the Jiangxi Central Station of the Irrigation Experiment. Daily meteorological conditions and 8 day NDVI at the Xiangtang station were shown in Figure 2.

Eddy Covariance Measurements and Meteorology Data at the Xiangtang Station
Up to date, the eddy-covariance technique has been the standard method in the ET and GPP measurements at scales of 10 2 -10 3 m [34]. To measure the fluxes of sensible heat, latent heat, and CO 2 , we have set up an open-path eddy covariance system (Li-COR 7500) at a height of 2.5 m since October 2016 at the Xiangtang station. The three-dimensional wind velocities and CO 2 /H 2 O concentrations were recorded by a datalogger (model CR3000, Campbell Scientific, Inc., Logan, UT, USA) at a frequency of 10 Hz.
Other microclimate variables were measured as 30 min averages, including net radiation (R n ) (CNR4, Kipp and Zonen, Delft, The Netherlands), air temperature (T a ), and relative humidity (RH) (HMP155A, Vaisala, Inc., Vantaa, Finland), wind speeds (U) and wind direction (WD) (model 03002, RM Young, Inc., Dundee, MI, USA), and precipitation (model 52202, RM Young, Inc.). The soil heat fluxes (model HFP01SC, Campbell Scientific, Inc.) were measured at 0.05 m below the soil surface. All data were collected by a CR3000 datalogger as 30 min averages. The data were selected from October 2016 to December 2018, except from 16 November 2017 to 15 March 2018, due to an instrument malfunction. The data were divided into 18 periods according to the growing stages of early rice and late rice, as shown in Table 1.  We used the EddyPro software v6.1.0 (www.licor.com/eddypro, accessed on 15 December 2019) to process the high-frequency raw data. After the 30-min scale fluxes were calculated using the high-frequency measurements in the air temperature, humidity, CO 2 , and three-dimensional wind velocity, the post-field data processing program was developed following FLUXNET's standard procedures [35]. The main data processing steps on the 30-min scale fluxes included spike removal [36], double coordinate rotation [37], frequency response corrections [38,39], and WPL corrections [40] for the air density fluctua-tions. The data quality check followed the method in Foken et al. [36], and the low-quality data were discarded in this study. Approximately 18.7% of the half-hour data points were removed as a result of a malfunction of the sensors, rain/fog influence, calibrations, and quality checks.
The net ecosystem exchange (NEE) measurements of carbon dioxide are then decomposed into GPP and R eco (Respiration of the ecosystem) using the nighttime-based partition method [41,42], in which the R eco = f(T a ) relationship is established using the nighttime NEE data. The daytime R eco is then estimated using the Reco = f(Ta) relationship and daytime T a , and therefore GPP during daytime is then estimated as the difference between daytime NEE and daytime R eco . The surface energy balance closure at the Xiangtang station was 71% on a daily scale. The heat storage correction on soil heat flux was calculated using the method of Heusinkveld et al. [43], and the energy balance correction on H and LE was performed using the Bowen-ratio method [44].
Data gaps were filled by three methods [45], including linear interpolation (gap < 2 h), look-up table method (gap within 2-12 h), and the multi-regression method (gap > 12 h). Missing net radiation, air temperature, relative humidity, and wind speed data were filled by data collected at the national standard meteorological station located near the Xiangtang station.

Regional Data
We used the China Meteorological Forcing Dataset (CMFD) [46,47] for regional climate forcing, which is at the spatial resolution of 0.1 • , and temporal resolution of 3 h. CMFD was produced by merging a variety of data sources, including CMA (China Meteorological Administration) weather station observation data, TRMM satellite precipitation analysis data (3B42), GEWEX-SRB downward shortwave radiation data, the Modern Era-Retrospective Analysis for Research and Applications (MERRA) (surface pressure), and GLDAS data (wind, air temperature, relative humidity, air pressure, and precipitation) [46].
The NDVI and albedo data used in this study were derived from the MODIS products MOD13Q1/MYD13Q1 and MCD43A3, which have a spatial resolution of 250 m and 1000 m, respectively. Due to the existence of cloud cover, the MODIS NDVI/Albedo values are often missing or underestimated. To remove the outliers and fill the missing gaps in the original time series, we used the technique of Harmonic Analysis of Time Series (HANTS) [48]. The LAI was estimated using the NDVI data. LAI is estimated from NDVI in this study due to its relatively high spatial resolution (250 m). First, vegetation coverage f c is estimated using the NDVI, i.e., f c = (NDVI max − NDVI)/(NDVI max − NDVI min ), where NDVI max and NDVI min are taken as 0.95 and 0.01, respectively. LAI is then inverted using the relationship between the LAI and vegetation coverage, i.e., f c = 1 − exp(−0.6 * LAI) [49,50]. Land cover data in the years 2000, 2005, 2010, and 2015 were provided by the National Earth System Science Data Center of China (http://auth.geodata.cn, accessed on 30 June 2020). All data were resampled to 250 m and a daily scale.

The PML-V2 Model
The PML-2 model is an updated version of the Penman-Monteith-Leuning (PML) model [22,23,51], which estimates ET using the PM equation with a Jarvis-Stewart type conductance model [52,53]. The Jarvis-Stewart type conductance model estimates the leaf conductance to CO 2 /vapor by multiplying the species-dependent maximum conductance values with scaling factors that reflect the constraints of radiation, leaf properties, and water supply conditions [50,54]. To make the conductance model more physically realistic, Gan et al. [30] replaced the semi-empirical Jarvis-Stewart model with a biophysical one [55] that reflects the control of photosynthesis, concentrations of CO 2 , and water vapor on leaf stomatal conductance.
The major components of the PML-V2 model, i.e., the estimation of GPP using the light-response curve, of ET using the PM model, and the coupling between GPP and ET using the biophysical conductance model, are introduced as follows. For more details, please refer to Gan et al. [30] and Zhang et al. [29].

The Estimation of GPP
In PML-V2, the photosynthesis rate of leaves, A g , estimated from radiation intensity (I, µmol m −2 s −1 ) and the concentration of CO 2 (C a , ppm) using a rectangular hyperbola function [10] are shown as follows: where β is the quantum efficiency, defined as the initial slope of the response curve of assimilation rate of leaves to light (µmol CO 2 (µmol PAR) −1 ), η is the carboxylation efficiency, defined as the initial slope of the response curve of assimilation rate of leaves to is the maximum photosynthetic rate obtained when both I and C a are saturated, which is temperature-dependent [56,57], shown as follows: where V m,25 is a parameter that has the same dimension with A m , and a and b are temperature coefficients, which are taken as 0.031 and 0.115, respectively [58]. The PML-V2 model estimates the canopy carbon assimilation rate A cg by integrating the A g of leaves to the canopy scales using LAI, shown as the following: where The ecosystem GPP is finally determined using A cg and a constraint function that reflects the effect of atmospheric dryness (D a , vapor pressure deficit, kPa) on GPP, where D max and D min are parameters.

The Estimation of ET
ET includes transpiration from leaf stomata (T), evaporation from the soil (E s ), and wet leaf surfaces (E i ) [59,60]. T is estimated using the Penman-Monteith model [22,23], shown as the following: where ε = ∆/γ(-), ∆ is the slope of the saturated vapor pressure to the air temperature (kPaK −1 ), γ is the psychrometric constant (kPa K −1 ), and A is the available energy (net radiation R n minus soil heat flux G, MJ m −2 d −1 ); ρ is the density of air (g m −3 ); c p is the specific heat of air at constant pressure (J g −1 K −1 ); G a is the aerodynamic conductance (m s −1 ), and G c is the canopy conductance (m s −1 ). R n was estimated as the difference between the incoming and outgoing radiation at the surface, i.e., where α is the albedo, R s is the incoming shortwave solar radiation, R li is the incoming longwave radiation, and R lo is the outgoing longwave radiation.
The soil evaporation (E s ) is estimated using the PT equation, shown as follows: where the fraction coefficient f varies from 0 to 1 when the soil surface is changing from completely dry to wet. In PML-V2, f is determined by the relative magnitude between the accumulated precipitation and soil equilibrium evaporation in the previous 32 days [29].
In this study, we added a parameter α s to account for the saturation property of the rice paddy field because a previous study showed that, on average, evaporation of wetland in this region is 1.26 times the equilibrium evaporation [61]: The evaporation from wet leaf surfaces (ET i ) is expressed as with where P wet is the reference threshold rainfall amount if the canopy is wet (mm d −1 ), f ER is the ratio of average evaporation rate over average precipitation intensity storms (unitless), f V is the fractional area covered by intercepting leaves (unitless), and P is the daily precipitation (mm d −1 ). S l and F 0 are two free parameters for E i . S l is the specific canopy rainfall storage capacity per unit leaf area, and F 0 is the specific ratio of average evaporation rate over average rainfall intensity during storms. For more details on the estimation of E i , please refer to the supplementary material of Zhang et al. [29].

Biophysical Conductance Model
The leaf conductance g s can be written as a function of A g , C a , and D a [55]: The PML-V2 model estimates the canopy conductance G c by integrating g s to the canopy scale, shown as follows:

Model Calibration and Validation
To evaluate the PML-V2 model, we calibrated the model parameters for each growing stage (P 3 , P 4 , . . . , P 9 , in Table 1) separately, in 2017, and tested the model performances in 2018, where parameters in each growing stage (P 11 , P 12 , . . . , P 17 , in Table 1) in 2018 were taken as those in the periods of P 3 , P 4 , . . . , P 9 , respectively, in 2017. The definitions and ranges of all parameters of the PML-V2 model are summarized in Table 2. For each period, we first calibrated the parameters of V m25 , β, η, D min , D max , and K Q by minimizing the RMSE in GPP predictions and then calibrated the rest of the parameters by minimizing the RMSE in ET predictions. In addition, to explore the seasonal variations in the model parameters that are related to the light responses, canopy conductances, and water-carbon coupling characteristics of the rice paddy ecosystem, we also calibrated the model in each period of early rice and late rice in both 2017 and 2018 (P 3 , P 4 , . . . , P 9 and P 11 , P 12 , . . . , P 17 ), separately. Table 2. Summary of parameters used in the PML-V2 model.

Parameters
Definition Ranges

V m25
The notional maximum catalytic capacity of Rubisco per unit leaf area at 25 • C 10-120 µmol m −2 s −1 β The initial slope of the response curve of assimilation rate of leaves to light (quantum efficiency) 0.01-0.15 µmol CO 2 (µmol PAR) −1 η The initial slope of the response curve of assimilation rate of leaves to CO 2 (carboxylation efficiency) 0.01-0.15 µmol m −2 s −1 (µmol m −2 s −1 ) −1

D min
The threshold below which there is no vapor pressure constraint 0.5-1.5 kPa D max The threshold above which there is no assimilation 3.5-6.5 kPa D 0 Water vapor pressure deficit of the air 0.5-2.0 kPa Extinction coefficient of available energy 0.7-0.9 S l Specific canopy rainfall storage capacity per unit leaf area 0.01-0.17

GPP and ET Fluxes of RICE Paddy Field at the Site-Level
As shown in Figure 3, the GPP of the rice paddy field at the Xiangtang station exhibited seasonal variations with two peaks, corresponding to the growth of the early rice and late rice throughout the year. For the early rice, the GPP peaks at the middle/late stages, i.e., 12.02 gCm −2 day −1 and 13.73 gCm −2 day −1 for 2017 and 2018, respectively ( Figure 3). On average, the GPP of early rice was 4.78 gCm −2 day −1 and 6.02 gCm −2 day −1 for 2017 and 2018, respectively (Table 3). In contrast, the average GPP of late rice, i.e., 8.70 gCm −2 day −1 and 7.91 gCm −2 day −1 for 2017 and 2018, respectively, were much larger than those for the early rice. In addition, compared to the early rice, GPP of the late rice increases more rapidly at the beginning of the growing stage ( Figure 3); therefore, GPP peaks at the middle stages for the late rice, and the peaks were larger than those for the early rice, i.e., 16.64 gCm −2 day −1 and 14.45 gCm −2 day −1 for 2017 and 2018, respectively. The faster growth rate in the beginning stage of the late rice than of early rice may be due to the fact that crops usually grow slower at lower temperatures. T a at the beginning stage (since 24 March) of the early rice was 16.7-21.8 • C lower than that at the beginning stage (since 15 July) of the late rice ( Figure 2).
Compared to the GPP, ET of the rice paddy has only one peak ( Figure 3). There appear to be two peaks in 2017; however, the low ET in late June and the first few days in July in 2017 were mainly due to cloudy and rainy conditions, in which the relatively small solar radiation (Figure 2) was probably the limiting factor for ET. As shown in Table 3, ET consumed most of the available energy because the Bowen ratio (the ratio between sensible and latent heat) of the paddy field was only about 0.1 for both the early rice and late rice. On average, ET for the late rice is generally larger than for early rice. However, the ecosystem water use efficiency (WUE), which is defined as the ratio between the GPP and ET, i.e., 1.86 gC/kgH 2 O and 1.59 gC/kgH 2 O for late rice, was still higher than those, i.e., 1.41 gC/kgH 2 O and 1.48 gC/kgH 2 O for early rice, in 2017 and 2018, respectively.

Seasonal Variations in the Key Parameters of the PML-V2 Model for the Early Rice and Late Rice Ecosystems
To reveal the seasonal variations in the key characteristics that govern the GPP, transpiration, and soil evaporation process, mainly including the light-response (β) and CO 2response (η) properties, the maximum catalytic capacity of Rubisco (V m25 ), the canopy conductance coefficient (m), and the soil evaporation coefficient (α s ), we calibrated the parameters of the PML-V2 model in every growing stage separately (Table 2) for both the early rice and the late rice ecosystems, as shown in Table 4.   The parameters that govern the GPP estimation differ for the early rice and late rice. The carboxylation efficiency η for the early rice increases as the crop grows, from less than 0.02 µmol m −2 s −1 (µmol m −2 s −1 ) −1 to nearly 0.15 µmol m −2 s −1 (µmol m −2 s −1 ) −1 at its grouting stage ( Figure 4) and then decreases. In addition, the maximum catalytic capacity of Rubisco, V m25, and the quantum efficiency β generally exhibited a similar change pattern as the GPP (Figure 4c) during the early rice period. Compared to the early rice case, the change patterns in parameters for the late rice periods did not typically coincide with the GPP. In addition, the inter-annual variations in parameters are larger for late rice. For example, the late rice in 2017 was less sensitive to light and CO 2 , with generally much smaller quantum efficiency β and the carboxylation efficiency η than those for late rice in 2018. The parameters that govern the transpiration estimation for the early rice and late rice showed similar averages and variability (Figure 4d,e). For example, the parameter m, which represents the sensitivity of the canopy conductance to the photosynthesis rate, is generally 20-25 on average for both the early rice and late rice. However, the parameter m showed smaller deviations from its mean value in the early rice period than in the late rice period. In addition, the parameter m generally decreases as GPP increases, indicating that the canopy conductance is more sensitive to GPP when GPP is smaller. The soil evaporation coefficient α s varies around 1.1 for both the early rice and late rice periods.

Model Evaluation at the Xiangtang Station
For the calibration period (2017), model performances were reasonably good, with predictions in the GPP and ET matching the dynamic characteristics and the peaks of measured GPP and ET ( Figure 5). The RMSE in GPP generally peaks in the stage of the resume growth before jointing; however, the RMSE is generally less than 50% of the mean GPP ( Figure 6). In addition, the mean bias in the GPP for the early rice is close to 0 on average and relatively small compared to the mean GPP during the growing period of early rice. Compared to GPP, the RMSE in ET predictions from the PML-V2 model is even smaller, about 10% of the mean ET. The overall good performances in the simulations of both GPP and ET indicate that the PML-V2 model is capable of predicting crop yield and its water use for rice paddies if the proper parameters are used.  To fully evaluate the PML-V2 model, we used the "calibration-validation" strategy and re-ran the model in 2018 using the parameters that were calibrated in 2017. The seasonal changes in GPP predictions generally coincided with the change patterns in the GPP measurements throughout the year in 2018; however, GPP was underestimated in the early rice stage (Figures 5a and 6). Compared to the "optimized" simulations, where parameters were optimized in every growing stage in 2017 and 2018 (Figure 6), RMSE in GPP increased~0.2-1.0 gCm −2 d −1 when the model was run at the early rice period, and the mean bias reached over −2 gCm −2 d −1 at the first growth stage of the early rice. In comparison, model simulations were much better for the late rice (Figures 5a and 6). Compared to the GPP estimation, ET estimation in 2018 showed more robust results. The predictions in ET and its three components were generally consistent with the optimized simulations. RMSE in ET predictions increased only~0.1-0.2 mm/day, and the mean bias was confined reasonably well.

Regional Estimation of GPP and ET Using the PML-V2 Model
We applied the PML-V2 model to estimate the regional (28-29 • N, 115.7-116.7 • E) distribution of GPP and ET in 2000, 2005, and 2015 at the spatial resolution of the NDVI data (0.0025 • ). The calibrated parameters, including those at different growing stages for the early rice and late rice at the Xiangtang Station, were used. Despite the relatively small extent of our study area (~10 4 km 2 ), some of the meteorological conditions exhibit large spatial variations even on an annual time scale, especially for the specific humidity, the wind speed, and the precipitation (Figures 7 and S1-S3). For example, annual precipitation ranges from 1400 mm to over 2200 mm across the region in 2005 ( Figure S2). In addition, meteorological conditions exhibit inter-annual variations (Figures 7 and S1-S3).   (Figure 8). We further demonstrated the spatial variability of GPP (ET and WUE) using the violin plot (Figure 9), which combines the boxplot of GPP (ET and WUE) with its probability density distribution. As shown in Figure 9, the total amount of GPP for the early rice and late rice also exhibits small inter-annual variations but large spatial variations within a single year. The probability density distribution of GPP was generally symmetric (Figure 9), indicating that the mean value of GPP is close to the median value of GPP. The total amount of GPP of the early rice was 709 ± 167 gC m −2 , 723 ± 174 gC m −2 , 681 ± 164 gC m −2 , and 707 ± 183 gC m −2 , and the total amount of GPP of the late rice was 605 ± 149 gC m −2 , 635 ± 167 gC m −2 , 638 ± 172 gC m −2 , and 746 ± 200 gC m −2 in 2000, 2005, 2010, and 2015, respectively. The GPP of the early rice showed less spatial variability than that of the late rice. Spatially, GPP was highly correlated (correlation coefficient r > 0.96) with the NDVI distribution for both the early rice and late rice periods ( Figure 10). In contrast, GPP for the early rice was weakly correlated with the meteorological conditions. For the late rice, GPP was negatively correlated with T a (r = −0.26).   The estimated annual ET of paddy rice (including early rice and late rice) was 988 ± 49 mm, 1075 ± 56 mm, 968 ± 46 mm, and 984 ± 48 mm in 2000, 2005, 2010, and 2015, respectively. No significant trend exists during the 2000-2015 period. Compared to the GPP, the estimated annual ET was also relatively large, but not always the largest, in the southeast of the region. For example, ET was the largest in the northern part of the region in 2005 due to the heavy rain in this area ( Figure S2). ET exhibits smaller spatial variations compared to that of GPP ( Figure 9). The probability density distribution of ET was also generally symmetric ( Figure 9). The total amount of ET of the early rice was 517 ± 19 mm, 551 ± 25 mm, 449 ± 15 mm, and 488 ± 22 mm, and the total amount of ET of the late rice was 471 ± 33 mm, 524 ± 34 mm, 518 ± 33 mm, and 496 ± 29 mm in 2000, 2005, 2010, and 2015, respectively. ET for the early rice was not strongly correlated with NDVI (r = 0.37). In contrast, ET was highly correlated with the solar radiation (r = 0.80) and air temperature (r = 0.77) for the early rice ( Figure 10). ET for the late rice was highly correlated with NDVI (r = 0.73), followed by wind speed (r = 0.54) and solar radiation (r = 0.48).
The median value of WUE is~1.40 gC/kgH 2 O for both the early rice and late rice. Despite the relatively small spatial variations in ET, the WUE exhibits larger variations than those of GPP and ET (Figure 9). WUE is generally higher in conditions with larger vegetation coverage, with correlation coefficients with NDVI equal 0.90 and 0.93 for the early and late rice, respectively. However, WUE is negatively correlated with R s , T a , u, and specific humidity.

Discussion
Rice is the major crop type in our study area. Compared to other crops, the WUE of rice is usually relatively low. For example, Wang et al. [62] reported that, among the staple crops, the maize cropland has the highest WUE with 2.48 ± 0.69 gC/kgH 2 O, followed by winter wheat (2.00 ± 0.39 gC/kgH 2 O), soybean (1.92 ± 0.52 gC/kgH 2 O), and paddy rice (1.88 ± 0.63 gC/kgH 2 O). Our results showed that WUE for rice at the Xiangtang station ranges from 1.41 gC/kgH 2 O to 1.86 gC/kgH 2 O for the early rice and late rice, which was consistent with the review of Wang et al. [62]. The relatively low WUE for rice is partly because the ET of rice is usually large due to the (nearly) constant saturation condition of the paddy fields. The Bowen ratio of the paddy field in our study was about only 0.1 for both the early rice and late rice ( Table 3), indicating that ET consumes most of the net radiation of the land surface, which is consistent with other studies, e.g., the Bowen ratio of the flooded rice paddy fields in the Philippines is also only 0.14 on average [63,64]. Among the total ET, transpiration takes up~50% on average in the growing stages of the early rice, and over 80% for the late rice ( Figure 11). Therefore, the ecosystem WUE (WUE e ) is close to the canopy WUE (WUE c ) for the late rice, whereas the WUE c is much larger than WUE e for the early rice, especially in the stages from jointing to yellow ripen ( Figure 11). The accuracy of the PML-V2 model in simulating GPP and ET of rice is reasonably good compared to other modeling studies. For example, Zhang et al. (2020) estimated the photosynthetic rate of rice using a light-response curve model and found that R 2 ranges from 0.77 to 0.96 at a daily scale when the multi-spectral reflectance data from an unmanned aerial vehicle were used. R 2 from the PML-V2 GPP simulations was also larger than 0.90 in our study case if parameters were optimized in each growing stage. R 2 was 0.72 for GPP during the validation period (2018) if the parameters in 2017 were used. As for ET estimation for paddy fields, Xu et al. (2017) proposed a Javis-Stewart type conductance model and successfully simulated rice ET using the conductance-based PM model at the Kunshan irrigation and drainage experiment station (31 • 15 15 N, 120 • 57 43 E), which is located to the northeast of the Xiangtang station (28 • 26 25 N, 116 • 00 01 E) and is also subject to a subtropical monsoon climate in the Taihu Lake region of China. After calibration, the model performed well with RMSE of 0.24 mm/d and 0.30 mm/d for 2014 and 2015, respectively, which is relatively small compared to the mean ET (3.52 mm/d) during the rice growth period (Xu et al., 2017). In comparison, the RMSE in ET predictions from the PML-V2 model is also less than 15% of the mean ET during the validation period (2018) if the parameters in 2017 were used.
Although similar modeling accuracy for GPP or ET has been achieved in previous studies, the advantage of the PML-V2 model is that the rates of photosynthesis and transpiration of vegetation are interconnected in the calculation process, which helps to correctly assess the WUE of an ecosystem and the canopy. However, we want to note that the reasonable performance of the model here relies on local calibration. In the paper of Zhang et al. (2019), parameters are fixed for all crops across the globe. Compared to their default values, e.g., the carboxylation efficiency (0.07 µmol m −2 s −1 (µmol m −2 s −1 ) −1 ), our study obtained similar values on average but with obvious seasonal variation. Other parameters, e.g., the stomatal conductance coefficient, m, varies around 20 in our study site, which is similar to that of the wetland (25.8) other than that of the crop (8.4) in Zhang et al. (2019).
We want to note that the PML-V2 model was evaluated using only one set of EC equipment at the Xiangtang site in our study area. However, because the study area is quite small, we assume that the parameters that were retrieved using the EC measurements at the Xiangtang site are representative of the characteristics of rice ecosystems in this area. Therefore, parameters of the PML-V2 model were assigned the same values for each rice grid in the study area. However, our results showed that GPP and ET of rice paddy fields exhibited substantial spatial variations (Figures 10 and 11) because meteorological and land surface conditions were different for different grids. Moreover, simple relations are expected to exist between the multi-year daily values of GPP, ET, and the meteorological and land surface variables. Compared to the meteorological variables which were at 0.01 • resolution, the NDVI data can be acquired at a much finer resolution (0.0025 • ), which was approximately equal to 250 m. Therefore, it is expected that the spatial distribution of estimated GPP was mainly correlated with the NDVI. The stepwise regression analysis showed that the NDVI can explain over 90% of the spatial variations in the multi-year daily values of GPP of the early rice. In addition, the slopes of NDVI in the regression equations for both the early rice and late rice were close to each other (14.14 v.s. 14.11). Compared to that of GPP, ET exhibited much smaller spatial variation due to the high evaporation rate of the saturated soil surface of paddy fields, which is comparable to the transpiration rate of the canopy; therefore, meteorological variables such as downward shortwave radiation, air temperature, and wind speed entered the regression equation. Besides GPP and ET, the stepwise regression also showed that, together with meteorological variables, the NDVI can explain over 90% of the spatial variations in the multi-year daily values of T/ET, and over 80% for those of WUE (Table 5).

Conclusions
Estimating GPP and ET in a water-carbon coupling way is important to correctly assess the water use efficiency of crops. In this paper, the usage of the PML-V2 watercarbon coupling model was evaluated in the lower reaches of the Poyang Lake plain. As a diagnostic model, the PML-V2 model does not require detailed crop management information; however, the model performed reasonably well in simulating both the GPP and ET and therefore the WUE because the key parameters of the model exhibited predictable seasonal variations, and the NDVI is capable of indicating crop growth status. Our study shows that the PML-V2 model is feasible in regional estimates of rice GPP and ET in the sub-tropical climate.
The vegetation information (e.g., NDVI) is the most important input for regional application because of its fine spatial resolution. However, the uncertainty of vegetation information will propagate into the estimation of GPP, canopy conductance, and therefore canopy transpiration. The calibrated parameters are prone to NDVI uncertainty. In addition, due to the limited temporal resolution of the MODIS NDVI product (16 days), it is difficult to retrieve the start and end dates of the growth stages of the rice. In conclusion, vegetation information at finer spatial and temporal resolution is the key to the application of the PML-V2 model.

Institutional Review Board Statement:
The study was conducted according to the guidelines of the Nanjing Institute of Geography and Limnology, Chinese Academy of Sciences, and approved by the Nanjing Institute of Geography and Limnology, Chinese Academy of Sciences.

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