Assessment of Potential Climate Change Effects on the Rice Yield and Water Footprint in the Nanliujiang Catchment, China

Abstract

The Nanliujiang catchment is one of major rice production bases of South China. Irrigation districts play an important role in rice production which requires a large quantity of water. There are potential risks on future climate change in response to rice production, agricultural irrigation water use and pollution control locally. The SWAT model was used to quantify the yield and water footprint (WF) of rice in this catchment. A combined method of automatic and manual sub-basin delineation was used for the model setup in this work to reflect the differences between irrigation districts in yield and water use of rice. We validated our simulations against observed leaf area index, biomass and yield of rice, evapotranspiration and runoff. The outputs of three GCMs (GFDL-ESM2M, IPSL-CM5A-LR and HadGEM2-ES) under three RCPs (RCP2.6, 4.5, 8.5) were fed to the SWAT model. The results showed that: (a) the SWAT model is an ideal tool to simulate rice development as well as hydrology; (b) there would be increases in rice yield ranged from +1.4 to +10.6% under climate projections of GFDL-ESM2M and IPSL-CM5A-LR but slight decreases ranged from −3.5 to −0.8% under that of HadGEM2-ES; (c) the yield and WFs of rice displayed clear differences in the catchment, with a characteristic that high in the south and low in the north, mainly due to the differences in climatic conditions, soil quality and fertilization amount; (d) there would be a decrease by 45.5% in blue WF with an increase by 88.1% in green WF, which could provide favorable conditions to enlarge irrigated areas and take technical measures for improving green water use efficiency of irrigation districts; (e) a clear rise in future grey WF would present enormous challenges for the protection of water resources and environmental pollution control in this catchment. So it should be to improved nutrient management strategies for the agricultural non-point source pollution control in irrigation districts, especially for the Hongchaojiang and Hepu irrigation districts.

1. Introduction

Over the past decades, the global climate has undergone significant changes with global warming [1]. Agriculture is the largest freshwater consumer, accounting for more than 70% of the world’s freshwater [2] and is extremely sensitive to climate change [3]. The change in the global climate will undoubtedly influence the future freshwater supply and food security by exerting pressure on agriculture and water resources, which presents a significant challenge to society [4]. It is necessary to research the impact of climate change on the utilization efficiency of agricultural water resources to mitigate the adverse impacts on agricultural production.

To assess water resource utilization during agricultural production processes, the water footprint (WF) concept has been proposed [5]. The calculation of a water footprint of a product can include the total amount of freshwater used along the supply chain of a product [6]. In the agricultural sector, WF of crop production has three components: (1) green WF: evaporation of water supplied from the rain during crop production; (2) blue WF: evaporation of the irrigation water supplied from groundwater or surface water sources; and (3) grey WF: volume of water required to dilute the pollution during the production process [7]. In recent years, various studies have focused mainly on the assessment of the water footprint in agricultural production on different scales (global, national, basin and district) from the perspective of production, climate change impact, or both [8,9,10].

Crop growth is a complex phenomenon that is closely related to water requirements and normally influenced by climate change. Rising temperature, elevated CO₂ and changing precipitation patterns are important aspects of the changing climate with pronounced impacts on crops in general [11]. (1) Effect of CO₂ increase on crop photosynthesis and respiration. Research has shown that the rates of crop photosynthesis and respiration increase along with the increment of the atmospheric CO₂ concentration as well as crop yield. Unlike C₄ crops (maize, sorghum, sugarcane), species of type C₃ (wheat, rice, soybean) have lower rates of respiration, and the efficiency of photosynthesis and its metabolic pathway react positively to high levels of CO₂ [8]. It had been reported that doubling the ambient CO₂ levels stimulates the biomass yield of C₃ plants by 40%, and data for C₄ plants indicate a stimulation of 11% [12]; (2) Effect of temperature change. Temperature controls the rate of crop metabolic processes that ultimately influence transpiration and grain production [13]. And high temperatures may increase the length of the growing period. It has been revealed that wild oat grown at a high temperature of 23/19 °C (day/night) complete their development faster than those grown at a normal temperature 20/16 °C [14]; (3) Effect of precipitation change. Crops show a particular sensitivity to the water supply. Decreased precipitation may lead to water stress, reducing the rain-fed crop output and increasing irrigation requirements in irrigated crops. However, the crop response to the interaction of CO₂, temperature and precipitation is complex. For instance, C₃ crops usually benefit from rising CO₂ levels physiologically, but rising temperatures can override the stimulating effects of CO₂ on the photosynthesis of C₃ crops [11]. Some studies have shown that high temperatures reduce or eliminate the high CO₂ growth enhancement [14].

Rice is one of the largest water consumers, utilizing 25% of the global agricultural freshwater supply [7]. Approximately 75% of rice production is obtained from irrigated rice [15]. Large irrigation projects are often constructed to meet the water demand in rice production. Additionally, rice is the third largest food crop planted, feeding more than 3 billion people worldwide [16]. According to the latest statistics reported by the USDA (U.S. Department of Agriculture), the world’s rice-sown area expanded to more than 160 million hectares by 2016, and the rice harvest exceeded 484 million metric tons [17]. Rice production has been chiefly concentrated in East, South and Southeast Asia, with a proportional planting area of 88.4% and yields of 89.3% yields worldwide. China has a rice sown area of 30.16 million hectares with a total rice output of approximately 145 million metric tons in 2016, making it the largest rice producer in the world [17]. Approximately 52% of the global rice harvested areas are significantly influenced by climate, including China, India and Indonesia with approximately 75% of the global rice production [18].

The tools used for modeling yield and WF of crops are expected to have process descriptions defining the crop growth and the stresses (water, temperature, and nutrient) associated with crop growth. Various crop models have been used to assess climate change in response to crop production and water use efficiency but with mixed results, such as the APSIM [19] model, the CropSyst [20] model, the DSSAT [21] model and the EPIC [22] model. However, these field-scale models were usually only used at specific sites, but limited to study regional scientific questions. The Soil and Water Assessment Tool (SWAT) is based on the inclusion of a crop growth module (EPIC) into a physically-based semi-distributed hydrologic model, which has been widely used to assess the impact of prospective climate change on the water resources and agricultural production of various basins across the world [23,24].

There are many works regarding yield and WF of upland crops production focused on dry-farming areas [25,26,27], but few quantitatively studies have been done on the yield and WF of rice production in South China, our understanding on climate change in response to rice production and water resources utilization in this area is still insufficient. In this study, the SWAT model is used to assess the yield and water footprint of rice in the Nanliujiang Catchment. The Catchment is an important double-cropping rice production region and makes an important contribution to grain consumed within South China. Irrigation districts play an important role in rice production, and cropping therein requires large amounts of water for irrigation, which is potentially at stake under future water stress. Furthermore, the provincial agricultural department plans to continue enlarging irrigated areas of rice paddy in the next 15 years. Therefore, the question arises whether future climate change may be potentially harmful to rice production, agricultural irrigation water management as well as pollution control locally. In addition, one needs a tool that has good adaptability to simulate crop production including fertilization and irrigation, and assess crop yield and water use under future climate change. Consequently, the aims of this study were (1) to evaluate performance of the SWAT model; (2) to evaluate the future climate of the Nanliujiang catchment; (3) to assess rice yield under current and future climate conditions; and (4) to investigate quantitatively the WF of rice production under current and future climate conditions.

2. Materials and Methods

2.1. Study Area

The Nanliujiang catchment (21°35′54″–22°52′32″N, 109°00′03″–110°23′12″E) is located in the south of Guangxi province, China. It borders the Darong Mountain on the North and the North Bay on the South. The 9569 km² catchment has a unique leaf-shaped physiographic characteristic with a drop of more than 1250 m, high-lying northeast, and low southwest (Figure 1). The Nanliujiang catchment belongs to the subtropics monsoon climate region, with an average annual temperature of 22 °C. The winter is mild with less rainfall, and the summer is hot and humid with abundant rainfall. The annual precipitation varies from 1596 mm in the northern mountains to 1804 mm on the southern coast. Having a proper soil and climate conditions, the catchment becomes one of the most important grain production bases in South China. The principal crops in this area are double-cropping rice. Early rice is usually planted in mid–March, until harvest in early July. Late rice is planted in mid–July and harvested around mid–November. The rainfall is heaving during the rainy season, about two–thirds of which falls during May–September; the maximum runoff usually occurs in July and August. However, the precipitation is less than the water requirement for rice growing in spring and autumn. Thus, irrigation is highly necessary for triggering rice growth and producing high yields. At present, the total cultivated land areas were expanded to 269.8 thousand hectares, approximately 70% of which were cultivated with rice (Figure 2a). There are 18 irrigating areas with a total area of 179.7 thousand hectares (Figure 2b) [28].

2.2. Description of the SWAT Model

The SWAT model is a physically semi-distributed hydrologic model, which has been widely used worldwide because of its powerful functionality for simulating water yield, sediment, water quality in basins [29,30], as well as predicting the impact of human activity and climate change on hydrology [31].

In the SWAT, the river basin is divided into subbasins, and each subbasin is further divided into several Hydrologic Response Units (HRUs). The hydrological component of SWAT calculates the daily water balance for each HRU and takes into account precipitation, evapotranspiration, soil moisture, surface runoff and percolation. The water balance equation is represented as

$$S{W}_t=S{W}_0+{\displaystyle \sum _{i=1}^T\left(P_{day}-R_{surf}-ET-W_{seep}-Q_{gw}\right)}$$

(1)

where SW is the soil water content; i is the time T (days) for the simulation period, P_day, R_surf, ET, W_seep and Q_gw are daily precipitation, surface runoff, evaporation, percolation and return flow (mm), respectively.

The plant growth module that is incorporated into SWAT is a biogeochemical-based EPIC model [32,33]. The daily heat unit accumulation is used to regulate plant growth. Plant growth and crop yield are constrained by temperature, water, atmospheric CO₂ concentration and fertilizer application.

Plant biomass production is modeled by simulating leaf area development, radiation interception and conversion of intercepted light into dry biomass with the radiation-use efficiency of a specified plant. The equations are expressed as follows:

$$bio={\displaystyle \sum _{i=1}^d\left\{RU{E}_i\cdot 0.5\cdot H_{day}\cdot \left[1-\exp \left(-k\cdot LAI\right)\right]\right\}}$$

(2)

where bio is the total dry biomass production in the growing season (kg·ha⁻¹); RUE_i is the radiation use efficiency on day i, kg·ha⁻¹·(MJ·m⁻²)⁻¹; H_day is the incident total solar (MJ·m⁻²); k is the light extinction coefficient; and LAI is the leaf area index.

The total leaf area index is calculated as follows:

$$LAI=\{\begin{array}{cc}{\displaystyle \sum _{i=1}^d\left(f{r}_{LA\mathrm{Im}x,i}-f{r}_{LA\mathrm{Im}x,i-1}\right)\cdot LA{I}_{mx}\cdot \left\{1-\exp \left[5 \left(LA{I}_{i-1}-LA{I}_{mx}\right)\right]\right\}}\hfill & f{r}_{PHU}\le f{r}_{PHU,sen}\hfill \\ LA{I}_{mx}\cdot \frac{1-f{r}_{PHU}}{1-f{r}_{PHU,sen}}\hfill & f{r}_{PHU}>f{r}_{PHU,sen}\hfill \end{array}$$

(3)

where fr_LAImx,i and fr_LAImx,i−1 are the fraction of the plant’s maximum leaf area index for day i and i − 1, respectively; LAI_mx is the maximum leaf area index of the plant; fr_PHU is the fraction of potential heat units accumulated for the plant on a given day in the growing season; fr_PHU,sen is the fraction of the growing season (PHU) at which senescence becomes the dominant growth process.

The crop yield is determined by the total crop biomass, and the harvest index and is calculated as follows:

$$Yield=bi{o}_g\cdot \left[\left(HI-H{I}_{\min }\right)\cdot \frac{\gamma }{\gamma +\exp \left(6.13-0.883\cdot \gamma \right)}+H{I}_{\min }\right]$$

(4)

where Yield is the crop yield (kg·ha⁻¹); bio_g is the aboveground biomass on the day of harvest (kg·ha⁻¹); HI is the potential harvest index on the day of harvest; HI_min is the minimum harvest index allowed for the plant; γ is the water deficiency factor.

The research shows that the change in the atmospheric CO₂ concentrations has a large influence on the radiation use efficiency [34]. The formulas expressing the relationship between the radiation-use efficiency and CO₂ concentrations have been incorporated into SWAT for climate change studies, and the equation is as follows:

$$RUE=\frac{100{\mathrm{CO}}_2}{{\mathrm{CO}}_2+\exp \left(r_1-r_2\cdot {\mathrm{CO}}_2\right)}$$

(5)

where CO₂ is the concentration of carbon dioxide in the atmosphere (ppm); r₁ and r₂ are shape coefficients; and RUE is the radiation use efficiency, kg·ha⁻¹⋅(MJ·m⁻²)⁻¹.

2.3. Data and SWAT Model Setup

The data required in this paper include the following: (a) a DEM map provided by the Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences; (b) a land use map (Figure 2a) and (c) a soil map obtained from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (d) daily meteorological data from 1995 to 2013 obtained from the China Meteorological Science Data Sharing Service Network, consisting of daily precipitation, maximum and minimum air temperature, relative humidity, sunshine duration and wind speed for 5 meteorological stations (Figure 1); (e) monthly pan evaporation data from 1995 to 2013 for 3 meteorological stations (Yulin, Bobai and Hepu) also obtained from the China Meteorological Science Data Sharing Service Network; (f) monthly measured streamflow data of three hydrometric stations (Hengjiang, Bobai and Changle) from 1995 to 2013; (g) the monthly outflow data of the 3 large reservoirs (Xiaojiang, Hongchaojiang and Wangshengjiang) during 1990–2013 provided by the Guangxi water conservancy and Electric Power Survey; (h) observed data of biomass and leaf area index for primary crops grown in the Hepu irrigation district from 2005 to 2007 provided by the Guangxi Academy of Agricultural Sciences, and rice yield data from 1999 to 2012 provided by the China Meteorological Science Data Sharing Service Network; (i) data for climate change scenarios from the output of three global climate models (

Table 1. Descriptions of the 3 regional climate change scenarios used in this study.

Content	Description
Model name	GFDL-ESM2M, HadGEM2-ES, IPSL-CM5A-LR
Scenario	RCP2.6, RCP4.5, RCP8.5
Spatial resolution	0.5° × 0.5°
Meteorological variables	Precipitation, maximum and minimum air temperature
Temporal scale	Daily data from 1 January 1960 to 31 December 2050

) provided by the Inter-Sectoral Impact Model Inter-comparison Project (ISI-MIP), which had been linearly interpolated and verified by statistical bias the probability distribution [35]. CO₂ was taken as a constant for each RCP (RCP2.6, CO₂ = 443, RCP4.5, CO₂ = 487 ppm, RCP8.5, CO₂ = 540 ppm, against CO₂ = 370 ppm currently) [8,24].

In plain agricultural regions with strong disturbance by human activity, the construction of ditch system projects profoundly alter the natural water flow paths and original sub-catchments [36,37]. Irrigation return water and runoff yield are mainly drained by the drainage canal system [38], making the irrigation district become an artificial subbasin [39,40]. Consequently, plain areas should not be divided into subbasins the watershed but the irrigation boundary controlled by drainage canal systems [41]. To show the real irrigation situation as closely as possible and reflect the differences between irrigation districts in yield and water use of rice, we used a combined method of automatic and manual delineation in this study. With the traditional watershed delineation method, the Nanliujiang catchment was divided into 26 subbasins according to the watershed (Figure 2c). Regarding the subbasins which contain both mountainous and plain, using manual delineation, 19 subbasins (1–2, 4–14, 16, 19–23) were redivided into 58 subbasins the boundaries of irrigation districts. Thus, a total number of 64 subbasins were identified for the catchment (Figure 2d), and 881 HRUs were defined therein. Depending on irrigation district location, irrigated agriculture was considered in subbasins 3, 6–7 and 24–59, and rain-fed cultivated land was considered in the others.

For each type of crop, the management practice was specified, such as planting, tillage, fertilization, irrigation and harvesting, among others. All operations were scheduled by the heat unit the cropping calendar. Fertilizers were applied with no effects of nitrogen and phosphorus stress through automatic fertilization application option in the SWAT model. The major sources of irrigation are reservoirs and river reaches (only irrigation district 4, 10 and 11 shown in Figure 2b). Irrigation water was extracted from the river reach or reservoir of a subbasin in response to the irrigation district location. Irrigation practices were accounted for in the simulations through automatic irrigation demand at the designated subbasins (it can be considered as sufficient irrigation without water stress), which aims to quantify the potential water requirement for rice growth under the influence of future climate change.

2.4. Crop Water Footprint

The water footprint (WF) is the total amount of fresh water needed to produce the goods and services related to a certain consumption pattern, which consists of green, blue, and grey WF [2]. The total WF of crop production can be calculated using the following equation:

$$WF=W{F}_{green}+W{F}_{blue}+W{F}_{grey}$$

(6)

where WF is the total WF of crop production (m³·kg⁻¹); and WF_green, WF_blue and WF_grey are the green, blue and grey WF (m³·kg⁻¹), respectively.

WF_green is calculated as follows:

$$W{F}_{green}=W_{green}/Y=10{\displaystyle \sum _{i=1}^{T}E{T}_{green}}/Y$$

(7)

where WF_green is the green WF (m³·kg⁻¹); W_green is the amount of green water use per unit area (m³·ha⁻¹); Y is the crop yield per unit area (kg·ha⁻¹); T is the length of the crop growing season (d); and ET_green is the daily green water evapotranspiration (mm), which is the minor factor between the daily crop evapotranspiration ET and the daily effective precipitation P_eff.

$$E{T}_{green}=\min \left(ET,P_{eff}\right)$$

(8)

where ET is the daily evapotranspiration (mm) and P_eff is the daily effective precipitation (mm).

Effective precipitation is the fraction of the total precipitation that is available to the crop and does not run off [42]. Based on this definition, P_eff is calculated by using the SWAT model as follows [6]:

$$P_{eff}=P-RO-DP$$

(9)

where P is rainfall, RO is surface runoff and DP is deep percolation.

W_blue is calculated as follows [2]:

$$W{F}_{blue}=W_{blue}/Y=10{\displaystyle \sum _{i=1}^{T}E{T}_{blue}}/Y$$

(10)

where WF_blue is the blue WF (m³·kg⁻¹); W_blue is the amount of blue water use per unit area (m³·ha⁻¹); and ET_blue is the daily blue water evapotranspiration (mm). In irrigated systems, ET_blue is the difference between ET and ET_green, and ET_blue = 0 in rain-fed systems [6,43].

$$E{T}_{blue}=\{\begin{array}{cc}ET-E{T}_{green}\hfill & in irrigated systems\hfill \\ 0\hfill & in rainfed systems\hfill \end{array}$$

(11)

WF_grey is defined as the freshwater required to assess the load of pollutants during crop growth, and it expresses the degradative water use [7].

$$W{F}_{grey}=a\cdot AR/\left(c_{\max }-c_{nat}\right)/Y$$

(12)

where WF_grey is the grey WF; AR is the product of the fertilizer application rate in the field per unit area (kg·ha⁻¹); and a is the leaching-runoff fraction. Nitrogen is the main pollutant considered for most water footprint studies, which is estimated due to its potential pollution effect via nitrate leaching into ground and surface waters, here a = 0.14 [43]; and c_max is the maximum acceptable concentration (10 mg·L⁻¹) and c_nat is the concentration in natural water (0.1 mg·L⁻¹) [44].

2.5. Simulation of Future Crop Yield and Water Footprint

To quantify the impact of climate change on the crop yield and WF of rice in the Nanliujiang catchment, the bias-corrected output of climate variables (precipitation, temperature, and CO₂) from GCMs were fed into the calibrated SWAT model, and thus rice yield and WF were calculated for future periods. Here, the reference period (2001–2010) was ten years, and therefore the future reference period of the same length (2041–2050) was used for improved comparability. It was assumed that there will be no change in land use for the project area in the future. Climate variables in the future, such as wind speed, relative humidity, and sunshine hours, were also assumed to be the same as the reference period.

3. Results and Discussion

3.1. Performance of the SWAT Model

3.1.1. LAI and Biomass

The crucial parameters that pertain to crop growth were calibrated on a daily time step to minimize the simulated and observed residuals for LAI and BIOM. The related parameters of early and late rice were displayed their initial and final values in

Table 2. Results from the calibration of SWAT rice growth parameters.

Parameter Name	Parameter Definition	Early Rice		Late Rice
		Initial Value	Final Value	Initial Value	Final Value
BIO_E	Biomass/energy ratio	22	24.87	22	28.71
HVSTI	Harvest index	0.5	0.54	0.5	0.55
FRGRW1	Fraction of the plant growing season corresponding to the first point on the optimal leaf area development curve (OLADC)	0.3	0.31	0.3	0.34
FRGRW2	Fraction of the plant growing season corresponding to the second point on the OLADC	0.7	0.51	0.7	0.38
LAIMX1	Fraction of the maximum leaf area index corresponding to the first point on the OLADC	0.01	0.28	0.01	0.34
LAIMX2	Fraction of the maximum leaf area index corresponding to the second point on the OLADC	0.95	0.99	0.95	0.58
BLAI	Maximum leaf area index	5	6.30	5	7.28
DLAI	Fraction of growing season when leaf area starts declining	0.8	0.79	0.8	0.76
EXT_COFF	Light extinction coefficient	0.35	0.57	0.35	0.56

. It can be seen that the values of BIO_E, FRGRW2, LAIMX2 and BLAI showed greater differences between early and late rice, which were caused mainly by rice varieties and climate variations in growing season. LAI and BIOM were adjusted observed measurements in the Hepu irrigation district from 2005–2007 (Figure 3). The performance indicators of calibration in this study were Nash–Sutcliffe (NS), root mean square error (PBIAS). Figure 3 and

Table 3. Model performance for the simulation of rice in the Nanliujiang catchment. PBIAS—Percent bias; NS—Nash-Sutcliffe coefficient.

Crop Type	Variable Name	Calibration (2005–2006)		Validation (2007)
		PBIAS	NS	PBIAS	NS
Early Rice	BIOM	−5.62	0.94	−7.43	0.97
	LAI	1.83	0.83	−3.82	0.87
Late Rice	BIOM	2.13	0.97	8.03	0.95
	LAI	−8.61	0.95	7.72	0.83

display that both early and late rice are well depicted, with PBIAS values ranging from −8.61 to 2.13 and NS values from 0.83 to 0.97. Similarly, validation indicated an acceptable model performance with PBIAS values ranging from −7.43 to 8.03 and NS from 0.83 to 0.97. It can be seen that the SWAT model can be used to simulate the development of rice well.

3.1.2. Discharge and Evapotranspiration

The first years (1995–1997) used to warm-up the model, and the data from 1998 to 2006 were used for calibration and those from 2006 to 2013 for validation. The 12 hydrological parameters were selected for 3 hydrological stations (Hengjiang, Bobai and Changle), and their optimal values for SWAT calibration are shown in

Table 4. Results from the calibration of SWAT hydrological parameters.

Parameter Name	Parameter Definition	Initial Range	Calibrated Value
			Hengjiang	Bobai	Changle
CN2	SCS runoff curve number	r (−0.5, 0.5)	0.10	−0.27	−0.24
ESCO	Soil evaporation compensation factor	v (0, 1)	0.21	0.24	0.31
EPCO	Plant uptake compensation factor	v (0, 1)	0.73	0.65	0.86
GWQMN	Threshold depth of water in the shallow aquifer required for return flow to occur	v (0, 5000)	1615.4	1248.5	1192.5
GW_REVAP	Groundwater “revap” coefficient	v (0.02, 2)	0.48	0.30	1.67
REVAPMN	Threshold depth of water in the shallow aquifer for “revap” to occur	v (0, 500)	105.67	334.0	164.39
SOL_AWC	Available water capacity of the soil layer	r (−0.5, 0.5)	0.28	0.27	0.28
SOL_K	Saturated hydraulic conductivity	r (−0.5, 0.5)	0.30	−0.14	0.17
ALPHA_BF	Base flow alpha factor	v (0, 1)	0.48	0.90	0.55
CANMX	Maximum canopy storage	r (−1, 1)	0.89	0.69	0.12
GW_DELAY	Groundwater delay	v (0, 500)	35.78	21.29	26.03
SLSUBBSN	Average slope length	r (−0.5, 0.5)	−0.26	−0.10	−0.25

Notes: r indicates the existing parameter value is multiplied by (1 + a given value), v indicates the default parameter is replaced by a given value.

. The results showed that the accuracy of the simulated value for the monthly streamflow at the Nanliujiang catchment in the calibration and verification period affectively met the standard (Figure 4). PBIAS of the simulated streamflow at the three stations ranged from 2.8 to 11.7 and from 3.5 to 13.0, respectively, and NS of the simulated streamflow at the three stations ranged from 0.85 to 0.94 and from 0.73 to 0.82, respectively (

Table 5. Model performance for the simulation of monthly discharge in the Nanliujiang catchment.

Station Name	Calibration (1998–2005)		Validation (2006–2013)
	PBIAS	NS	PBIAS	NS
Hengjiang	11.72	0.85	13.02	0.73
Bobai	3.22	0.90	3.54	0.87
Changle	2.81	0.94	4.33	0.82

).

A direct relationship has been between crop development and evapotranspiration [45], and thus the evapotranspiration is also necessary for calibration in this work. However, the data of the actual evapotranspiration is hard to be measured or recorded, and therefore the measured monthly pan evaporation data from 2000 to 2011 for 3 stations (Yulin, Bobai and Hepu) may be a good source of information. Here, the values of the actual evapotranspiration were estimated by using a simple conversion method. Xie et al., 2009 [46] examined the ratio between the actual evapotranspiration and pan evaporation in the Dongjiang catchment, which is near to the Nanliujiang catchment and shows similar climatic characteristics. The actual evapotranspiration was estimated by multiplying the measured pan evaporation with the ratio. The calibration results showed an acceptable performance of the model with PBIAS values from −4.2 to 10.4 and NS values from 0.70 to 0.85. Similarly, validation indicated an acceptable model performance with PBIAS values ranging from −8.42 to 13.2 and NS from 0.59 to 0.79 (

Table 6. Model performance for the simulation of monthly evapotranspiration in the Nanliujiang catchment.

Station Name	Calibration (2000–2005)		Validation (2006–2011)
	PBIAS	NS	PBIAS	NS
Yulin	−4.21	0.70	−6.53	0.59
Bobai	10.4	0.74	13.2	0.61
Hepu	−7.13	0.85	−8.42	0.79

).

3.1.3. Rice Yields

Root mean squared errors (RMSE) and mean relative error (MRE) were used to evaluate the SWAT performance. The parameters were calibrated measured yield data from two stations (Yulin and Hepu) for the period from 1999–2012.

Table 7. Model performance for the simulation of rice yield.

Station Name	Crop Type	Calibration (1999–2006)			Validation (2007–2012)
		PBIAS	RMSE/ton	MRE/%	PBIAS	RMSE/ton	MRE/%
Yulin	Early Rice	−13.13	0.83	13.11	−8.55	0.59	6.11
	Late Rice	7.90	0.63	11.42	1.49	0.62	8.98
Hepu	Early Rice	0.31	0.24	7.15	4.58	0.29	8.61
	Late Rice	−6.63	0.56	6.83	2.01	0.49	5.02

showed acceptable performance by the model, with PBIAS values ranging from −13.13 to 7.90 and RMRE values from 0.24 to 0.83, and the MRE between the observed and simulated yield was less than 15%. All the results were considered acceptable for calibration of the model. Figure 4 showed the simulated and measured rice yields at Hepu station from 1999 to 2012.

3.2. Crop Yield and Water Footprint in the Reference Period

The rice yield in the Nanliujiang catchment was estimated using the validated SWAT model for the reference period (2001–2010), as shown in Figure 5a. The estimated total, green, blue and grey WF of rice in the 2000s calculated by Equations (6)–(12) is shown in Figure 5b–e respectively. Statistics (mean) of the yield and WF values are also shown in the figures. A combination of water use, soil and fertilization amount determined the spatial distribution characteristics of rice production in the catchment. Figure 5a shows that the yields of irrigated rice were higher than those of rain-fed rice. Since the rice in the rain-fed area was grown under rain-fed conditions during the growing period, there was a rapid decline in yield relative to that grown in irrigation districts under the same climate conditions. Rice yields varied from 4761 to 5264 kg·ha⁻¹ in the irrigated area and from 3035 to 3378 kg·ha⁻¹ in the rain-fed area. In the irrigated subbasins, the rice yields in the southern catchment were usually larger than that in northern catchment, and the rain-fed subbasins shared a similar characteristic. Hepu and Hongchaojiang irrigation districts (subbasin 24, 48–59) showed the maximum yield per acre in all irrigation districts, vary from 5174 to 5264 kg·ha⁻¹. The Maolin irrigation district (subbasin 31) showed the minimum yield per acre at 4761 kg·ha⁻¹.

The average total WF of rice production in the study area was 1064 m³·t⁻¹ (Figure 5b), which was higher than the national average (971 m³·t⁻¹) but below the global average (1450 m³·t⁻¹) [47]. The largest share of the total WF corresponded to the green WF, accounting for 44.8%, followed by blue and grey WFs, which accounted for 29.1% and 26.1%, respectively (Figure 5d,e). The Hongchaojiang irrigation district (subbasin 56 and 59) showed the maximum WF in all irrigation districts of 1049 m³·t⁻¹, with the green WF representing 179 m³·t⁻¹, blue 591 m³·t⁻¹, and grey 279 m³·t⁻¹; the Longmen irrigation district (subbasin 27) showed the minimum total WF, consisting of a green WF of 200 m³·t⁻¹, blue of 490 m³·t⁻¹, and grey of 213 m³·t⁻¹, resulting in a total WF of 903 m³·t⁻¹. Subbasin 27 showed the minimum WF of 903 m³·t⁻¹, and subbasin 9 showed the maximum WF of 1248 m³·t⁻¹. The values of the WF of rain-fed subbasins with low rice production were higher than those of irrigated subbasins, confirming that a greater reduction in yield can result in large water footprints [6].

3.3. Future Climate Change

Since the south long-narrow catchment experiences a spatial and temporal variation in precipitation and temperature, the projection of future climate differs in different stations, which can be observed in the results provided by the different GCMs and RCPs. Future average precipitation and maximum and minimum temperature in the 2040s are shown in

Table 8. Average values of climate variables during the growth season for the 2000s and 2040s under 9 climate change scenarios.

Period	Scenarios		Precipitation (mm)		Temperature (°C)
					Maximum		Minimum
			Yulin	Hepu	Yulin	Hepu	Yulin	Hepu
2000s	RE		1390	1705	30.6	29.8	22.8	23.5
2040s	GFDL-ESM2M	RCP2.6	1553	2043	29.6	30.5	22.3	24.0
		RCP4.5	1647	1953	29.4	30.4	22.5	24.2
		RCP8.5	1377	1655	30.6	31.6	23.1	24.9
	HadGEM2-ES	RCP2.6	1364	1605	31.3	32.2	23.3	25.3
		RCP4.5	1423	1785	31.2	32.0	23.5	25.4
		RCP8.5	1450	1771	31.5	32.3	23.8	25.6
	IPSL-CM5A-LR	RCP2.6	1420	1706	30.2	31.2	22.9	24.6
		RCP4.5	1451	1871	30.4	31.3	23.2	25.0
		RCP8.5	1413	1731	31.0	32.0	23.7	25.5

during the growth season (March to November), under reference period (RE) and climate change scenarios. RCP8.5 of GFDL-ESM2M and RCP2.6 of HadGEM2-ES show a slight decrease in precipitation, whereas the others show an increase from 1390 mm currently to 1647 mm under RCP4.5 of GFDL-ESM2M at Yulin station, and from 1705 mm currently to 2043 mm under RCP2.6 of GFDL-ESM2M at Hepu station; RCP2.6 and 4.5 of GFDL-ESM2M show a decrease in temperature, whereas the others show an increase from 30.6 °C to 31.5 °C at a maximum temperature and from 22.8 °C to 23.8 °C at a minimum temperature under RCP8.5 of HadGEM2-ES at Yulin station, and from 29.8 °C currently to 32.3 °C at a maximum temperature and 23.5 °C currently to 25.6 °C at a minimum temperature under RCP8.5 of HadGEM2-ES at Hepu station.

3.4. Rice Yield and Water Footprint in the Future Climate

3.4.1. Rice Yield in the Future Climate

Rice is a C₃ plant, and therefore an increase in CO₂ can, in theory, increase its yield; however, yield is controlled by a combined effect of precipitation, temperature and CO₂. A higher temperature can shorten the period of rice growth (according to the heat unit approach), and temperature stress can reduce the accumulation rate of dry materials, both of which may lead to low production. Figure 6 shows the potential rice yield in the 2040s under 9 climate change scenarios. On average, rice yield would increase by +3.0% in the future relative to the 2000s, ranging from −3.5% (HadGEM2-ES, RCP8.5) to +10.5% (GFDL-ESM2M, RCP4.5). Both GFDL-ESM2M and IPSL-CM5A-LR provide an increasing yield by different degrees under all RCPs, but HadGEM2-ES instead depict a decreasing yield under all RCPs (Figure 6). GFDL-ESM2M shows an increased yield up to +10.9% for irrigated rice and +10.3% for rain-fed rice under RCP4.5 in the 2040s relative to the 2000s, due to increasing precipitation and CO₂. HadGEM2-ES provide the lowest rice yield due to the highest temperatures compared with the other models under all RCPs (Table 8), because rice yield will be largely impacted by temperature stress as well as a shortened growing period.

The spatial distributions of the average rice yields of the 3 GCMs under 3 RCPs in the 2040s were shown in Figure 7. It can be seen that the spatially clustered areas with high rice production are distributed in the northern (subbasin 3, 6, 7, 25–39), middle (subbasin 40–47) and southern (subbasin 24, 48–59) study area, which are the important irrigation districts and grain production base in South China. This finding can be attributed to strong factors such as a smooth topography, favorable soil quality and proper irrigation facilities. The average rice yield in the irrigation districts will be increased by 239 kg·ha⁻¹ under RCP2.6, 380 kg·ha⁻¹ under RCP4.5 and 199 kg·ha⁻¹ under RCP8.5 relative to the 2000s, respectively. Subbasins with low rice yield are distributed in mountainous and hilly areas with rain-fed farming due to a shortage of water conservancy facilities. The mountainous subbasins with larger yields will be subbasin 14–16, 18–21 and 64, which profit from advantages in climate conditions. The average rice yield will be increased by 138 kg·ha⁻¹ under RCP2.6, 223 kg·ha⁻¹ under RCP4.5 and 105 kg·ha⁻¹ under RCP8.5 relative to the 2000s, respectively. However, their rice yields will be no larger than 3400 kg·ha⁻¹. These mountainous subbasins show a rugged topography with poor soils, mountainous cultivation could not only destroy mountain vegetation but also lead to soil erosion and environmental pollution. Therefore, reasonable adjustment and optimization of the agricultural structure is of great necessity in this area in the future. Mountain areas may be more adaptable to planting economic forests that can generate better economic returns, such as longan, litchi, mango and citrus.

3.4.2. Water Footprint in the Future Climate

Figure 8a shows that the average total WF of irrigated rice in the 2040s will decrease by different degrees under all scenarios. On average, the total WF will be reduced by 2.3% under RCP2.6, 4.9% under RCP4.5 and 5.4% under RCP8.5 relative to the 2000s, respectively. A projected rise in rice yield along with a decrease in blue water may be the main cause of this phenomenon. Due to the rise in precipitation during the growing season, rice has more green water to meet the increasing evapotranspiration demand triggered by climate change, leading to a substantial increase in green water demand and a reduction in blue water use. Hence, blue WF shows a significant decrease by 43.0% under RCP2.6, 48.3% under RCP4.5 and 45.2% under RCP8.5, respectively. Accordingly, green WF shows an increase from 207 m³·t⁻¹ to 394 m³·t⁻¹ under RCP2.6, 395 m³·t⁻¹ under RCP4.5 and 382 m³·t⁻¹ under RCP8.5, respectively. The amount of nitrogen fertilizer application will significantly increase with rice yield growth in the future, whereas the gray WF is likely to slightly increase by 4.0% under RCP2.6, 4.1% under RCP4.5 and 1.0% under RCP8.5, respectively, because rice yield improvements have helped to reduce obvious growth in the gray WF. Figure 8b shows that the mean total WF of rain-fed rice under 9 scenarios in the 2040s showed no obvious upward or downward trend. A slight decline in green WF and rise in gray WF may be the main causes. The green WF shows a small decrease by 4.0% under RCP2.6, 6.7% under RCP4.5 and 7.0% under RCP8.5, respectively. The amount of fertilizer applied to rain-fed rice will be important for growth along with rice yield growth in the 2040s, and the gray WF is likely to increase under all scenarios by 13.9% under RCP2.6, 17.1% under RCP4.5 and 10.2% under RCP8.5, respectively.

The spatial distributions of the average total, green, blue and gray WF of the 3 GCMs under 3 RCPs in the 2040s are shown in Figure 9. The WFs of rice production exhibit cleat differences among subbasins under 3 RCPs, mainly due to the differences in rice yields and water use levels. The rice yields and water use differ between subbasins and irrigation districts because of their variances in climatic conditions, soil quality and fertilizing amount caused by climatic conditions. Although the values of total WFs are different under the 3 RCPs as well as the green, blue and gray WFs, their distributions under different RCPs are similar (Figure 9a vs. Figure 9e vs. Figure 9i for total WF; Figure 9b vs. Figure 9f vs. Figure 9j for green WF; Figure 9c vs. Figure 9g vs. Figure 9k for blue WF; Figure 9d vs. Figure 9h vs. Figure 9l for gray WF). For instance, Figure 9a,e,i show that the subbasins in the northern, middle and southern catchment have a lower green, gray and total WF, representing the subbasins located in plains planted with irrigated rice. The irrigation districts located in the southern catchment usually have higher WFs than other irrigation districts due to the high precipitation and temperature. Hongchaojiang irrigation district (subbasin 56, 59) and Longmen irrigation district (subbasin 27) show the maximum and minimum total WF among all the irrigation districts, respectively. The spatial distribution of green WF is consistent with the distribution of precipitation. Due to abundant rainfall, the proportion of green WF in the Hepu (Subbasin 24, 48–55 and 57–58) and Hongchaojiang irrigation districts is relatively high. Farmers should take technical measures to increase green water use and reduce irrigation water use as more as possible. The southern irrigation districts also have the largest gray WF in the catchment, which increase by 10.8% under RCP2.6, 13.2% under RCP4.5 and 8.3% under RCP8.5, respectively. It will increase the pressure of water resource protection and water pollution control in the downstream area. So the fertilizers application should be regulate scientifically in the future so that grey WF can be kept at a proper level.

4. Conclusions

In this study, we used the SWAT model to quantify the yield and WFs of double cropping rice within the Nanliujiang catchment, China. We validated the simulations against observed leaf area index, biomass, runoff, evapotranspiration and rice yield in sequence. The impact of potential climate change on the yield and WF of rice was analyzed in the 2040s under 3 GCMs with 3 RCPs. The results showed that: (a) the SWAT model is an ideal tool to simulate rice development and WF under given and future climate conditions; (b) yield simulation showed significant variations under different future climate scenarios. It showed increased rice yields ranged from +1.4 to +10.6% under the GFDL-ESM2M and IPSL-CM5A-LR models but slightly decreased yields ranged from −3.5 to −0.8% under the HadGEM2-ES model; (c) both the yield and WF of rice exhibited obvious differences among subbasins, with characteristically high values in the south and low values in the north, mainly due to differences in climatic conditions, soil quality and fertilization amounts; (d) blue WF would significantly decrease by 45.5% and green WF would show a great increase by 88.1%, meaning green water would contribute the largest part in the WF of rice production in the catchment. These findings would be good news for local water resource development and utilization, and it would provide favorable conditions to expand the area under irrigation and take technical measures for improving green water use efficiency in irrigation districts; (e) the gray WF would increase under potential climate change, suggesting an increased requirement for fertilizer application in the future. The finding suggested that it was necessary to improved fertilization management strategies for the agricultural non-point source pollution control in irrigation districts, especially for the Hongchaojiang and Hepu irrigation districts. This work will also delineate some helpful information for the management and protection of water resources as well as irrigation planning.