*Article* **Water Dynamics and Hydraulic Functions in Sandy Soils: Limitations to Sugarcane Cultivation in Southern Brazil**

**Jessica Lima Viana 1,\*, Jorge Luiz Moretti de Souza <sup>2</sup> , André Carlos Auler <sup>2</sup> , Ricardo Augusto de Oliveira <sup>3</sup> , Renã Moreira Araújo <sup>3</sup> , Aaron Kinyu Hoshide 1,4, Daniel Carneiro de Abreu 1,5 and Wininton Mendes da Silva <sup>6</sup>**


**Abstract:** Crop cultivation on sandy soils is susceptible to water stress. Therefore, we determined the physical-hydric attributes of a Latossolo Vermelho distrófico (Oxisol) in northwestern Paraná state, Brazil. Soil samples were collected at depth ranges of 0 to 0.2 m, 0.2 to 0.4 m, and 0.4 to 0.6 m. We measured clay, silt, sand, fine and coarse sand contents, soil particle density, soil bulk density, total porosity, microporosity, and macroporosity. We also measured soil characteristics such as saturated and unsaturated soil hydraulic conductivities, pore distribution, water retention, available water capacity, and easily available water. We also estimated soil moisture, matric potential at field capacity, and time at field capacity. Validation of associations among these soil physical-hydric attributes was performed using principal component analysis. For the sandy soils analyzed, the distributions of coarse and fine sand fractions were measured for better evaluation of the soil's physical and hydric attributes. Higher coarse sand contents increased soil hydraulic conductivities, maximum pore diameter, and macroporosity while reducing microporosity. Fine sand content reduced conductivity and increased soil water retention in subsurface layers. Simulated sugarcane yield increased with soil water storage. These results support improving crop simulation modeling of sugarcane to support sustainable intensification in regions with sandy soils.

**Keywords:** available water capacity; hydraulic conductivity; pore distribution curve; retention curve; soil texture

### **1. Introduction**

Physical-hydric attributes are used to understand and evaluate soil variability and quality, aeration, hydraulic conductivity, water redistribution, storage capacity, water availability to plants, and root growth [1–3]. Soil density, texture, total porosity, and saturated soil hydraulic conductivity are physical-hydric attributes used as primary indicators to monitor soil quality [4]. Sandy soils are associated with high hydraulic conductivity, and higher proportions of the coarse sand fraction are associated with reduced water availability [5,6].

In addition to intrinsic factors, soil use and management interfere with soil water retention [7] and consequently water availability for plants. A previous study on soil developed from sandstone in the 0 to 0.4 m layer, found 48.6 mm and 58.7 mm of available water capacity (AWC) in soil under crop-livestock integration and pineapple, respectively [6]. This difference of 10.1 mm was influenced by land use and management. Due to the intense

**Citation:** Viana, J.L.; de Souza, J.L.M.; Auler, A.C.; de Oliveira, R.A.; Araújo, R.M.; Hoshide, A.K.; de Abreu, D.C.; da Silva, W.M. Water Dynamics and Hydraulic Functions in Sandy Soils: Limitations to Sugarcane Cultivation in Southern Brazil. *Sustainability* **2023**, *15*, 7456. https://doi.org/10.3390/ su15097456

Academic Editor: Jan Hopmans

Received: 31 January 2023 Revised: 22 April 2023 Accepted: 25 April 2023 Published: 1 May 2023

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

management of the sugarcane crop, knowledge of the physical-hydric attributes of the soil is important, especially for crops established on sandy soils. The effect of sugarcane cultivation can increase soil bulk density, reaching critical values of 1760 and 1770 kg/m<sup>3</sup> in layers 0 to 0.2 m and 0.2 to 0.4 m, respectively, in medium texture Oxisol [8].

In sandy soils, there are differences in soil water retention due to the particle size composition of the sand fraction [9]. High proportions of medium and fine sand in sandy soils favor the formation of a capillary distribution network with smaller diameter pores, allowing water retention between soil particles. This matrix provides less soil moisture loss from percolation and, consequently, greater water storage due to the high percentages of finer particles (between 40% and 75% of total sand), combined with low clay content (i.e., particle diameters < 0.02 mm). Although typically less of a contributing factor, clay content can also contribute to increased water retention in sandy soils [9].

The Experimental Station of the Sugarcane Genetic Improvement at the Federal University of Paraná, in collaboration with the Inter-University Network for the Development of Sugarcane Energy, annually conducts numerous studies aimed at genetic improvement and productivity of sugarcane [10,11]. The results obtained in the experiments developed at this Experimental Station are a reference for numerous agricultural planning activities in the northwest region of Paraná state in southern Brazil. However, no detailed study of the physical-hydric attributes has yet been conducted for soils on the Experimental Station that represent a large number of Oxisols in the region. Obtaining information about the dynamics of water in the soil, especially because it is a region dominated by soils developed in the Caiuá Sandstone, may help in decision making regarding the selection of sugarcane varieties better adapted to water stress conditions. The survey of physical and hydric attributes can also subsidize water management in irrigated areas in sandy soils in northwestern Paraná. In this context, the goals and objectives of our research were to (1) determine the main physical-hydric attributes of the Latossolo Vermelho distrófico (Oxisol) cultivated with sugarcane in the northwestern region of Paraná state in Brazil, and (2) evaluate this soil's physical qualities for supporting the sugarcane crop.

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

#### *2.1. Study Site Characteristics*

This work was developed in the Laboratory of Modeling of Agricultural Systems (LAMOSA) at the Setor de Ciências Agrárias (SCA) at the Federal University of Paraná (UFPR). The analyses were performed using soil samples taken from the Experimental Station of the Sugarcane Genetic Improvement Program (PMGCA), belonging to UFPR and to the Inter-University Network for the Development of Sugarcane Energy (RIDESA). The soil of the Paranavaí Experiment Station is classified as Latossolo Vermelho distrófico (Oxisol) [12]. The PMGCA is located at 22◦580 South latitude, 52◦280 West longitude and is at 470 m of average altitude in the municipality of Paranavaí, located in the northwestern region of the Paraná state, Brazil. The climate of the region, according to the Köppen classification is Cfa (subtropical climate) and has an average annual air temperature between 22.1 and 23 ◦C, with average annual precipitation between 1400 and 1600 mm [13].

The geomorphology of northwestern Paraná is characterized by the Caiuá sandstone, with an area of 3.2 million hectares, representing 16% of the state area [14]. The Caiuá formation consists of deposits of eolic and fluvial environments, represented by fine to medium sandstones that are purplish with large cross-stratification [15]. According to Costa et al., 2020, the low natural fertility is due to the high sand content associated with the predominance of minerals of low activity in the clay fraction [16]. Latossolo soils associated with the Caiuá sandstone can also be highly susceptible to erosion as well as vertical subsurface soil water movement [17]. According to Pereira 2016, the Caiuá Aquifer is an important water source for the northwest region of the state of Paraná, with an outcrop area of approximately 30,000 km<sup>2</sup> . This aquifer is characterized by being free and porous, with bicarbonated calcic to calcium-magnesian waters. The aquifer is composed of aeolian

sandstones from the Caiuá Group, which were deposited during the Upper Cretaceous over a paleo-depression formed by the basaltic flows of the Serra Geral Formation [18].

#### *2.2. Soil Attributes*

In order to determine soil physical-hydric attributes, we collected soil samples at ten representative points in the experimental area during July 2019. Interval soil sampling was conducted for the soil layers at depths of 0 to 0.2, 0.2 to 0.4, and 0.4 to 0.6 m. Soil analyses were performed using the Soil Physics Laboratories of the Soil and Agricultural Engineering Department at UFPR and in the Rural Development Institute (IDR) in Londrina, Paraná state in Brazil.

The deformed samples were collected with a flat bottom scoop, dried in an oven at 40 ◦C, sieved through a 2 mm mesh, and stored in plastic bags. After that, the clay, silt, and sand contents were gravimetrically determined using the pipette method [19]. Separation of the total sand fraction into coarse sand (0.2 to 2 mm) and fine sand (0.05 to 0.2 mm) was also performed. These grain fractions were then weighed (g/kg), according to the methodology developed by Embrapa [19]. Soil particle density (ρP; kg/m<sup>3</sup> ) was determined by the modified volumetric balloon method [20] using the formula:

$$\rho\_{\rm P} = \frac{\mathbf{M}\_{\rm vs} - \mathbf{M}\_{\rm V}}{50 - \left(\frac{\mathbf{M}\_{\rm vs} - \mathbf{M}\_{\rm vs}}{\rho\_{\rm a}}\right)} \tag{1}$$

where ρ<sup>P</sup> is the particle density of soil (kg/m<sup>3</sup> ), Mvs is the mass of the volumetric flask containing soil (kg), M<sup>v</sup> is the mass of the volumetric flask (kg), Mvsa is the mass of the volumetric flask containing soil plus alcohol (kg), ρ<sup>a</sup> is the density of the alcohol (kg/m<sup>3</sup> ), determined by the relationship between M<sup>a</sup> and Vv, where M<sup>a</sup> is the mass of alcohol (kg), and V<sup>v</sup> is the volume of the flask (m<sup>3</sup> ).

#### *2.3. Saturated Soil Hydraulic Conductivity*

The un-deformed soil samples were collected in duplicate at each sampling point with the aid of volumetric rings (4.69 cm internal diameter and 3.51 cm height) and an Uhland sampler. After adequate preparation, the samples were saturated using capillary rise. After saturation, the saturated soil hydraulic conductivity (K<sup>S</sup> measured as cm/hour) was determined by the Constant Load Permeameter method [21]:

$$\mathbf{K}\_{\rm S} = \frac{\mathbf{V} \cdot \mathbf{H}}{(\mathbf{H} + \mathbf{h}) \cdot \mathbf{A} \cdot \mathbf{t}} \tag{2}$$

where K<sup>S</sup> is the saturated hydraulic conductivity (cm/hour), V is the volume of percolated water (cm<sup>3</sup> ), H is the height of the cylindrical ring with soil (cm), h is the water sheet over the sample (cm), A is the cross-sectional area of the sample (cm<sup>2</sup> ), and t is the percolation time (hours).

Total porosity (TP), measured in cubic meters (m<sup>3</sup> ) of air space per volume of soil (m<sup>3</sup> ), was considered equal to the volumetric soil moisture at saturation (θs), which was also measured in m3/m<sup>3</sup> of soil. The microporosity measured in m3/m<sup>3</sup> of soil corresponded to the volumetric moisture of the soil sample submitted to a tension of −6 kilopascals (kPa). The macroporosity also measured in m3/m<sup>3</sup> of soil was obtained by taking the difference between the value for TP and the value for microporosity [19].

#### *2.4. Soil Water Retention*

The soil samples were submitted to tensions of −6 and −10 kPa in a tension table and −33, −100, −500, and −1500 kPa in a Richards chamber. After each applied matric potential, the samples were dried at 105◦C for 48 h for the determination of soil bulk density and volumetric moisture (θ), according to the methodologies developed at Embrapa [19]. The "θ vs. ψm" points and the moisture at saturation were used to determine the soil water retention curve (SWRC), fitted with the model of van Genuchten 1980 [22] using the Burdine

1953 restriction [23]. The adjustment of the parameters of the van Genuchten model and estimation of the volumetric soil moisture (θˆ) was performed with the *nls* function of the open-course statistical software program R version 4.0.0 [24]:

$$\Theta = \frac{\theta(\psi\_{\rm m}) - \theta \mathbf{r}}{\theta \mathbf{s} - \theta \mathbf{r}} = \frac{1}{\left[1 + |\alpha \cdot \psi\_{\rm m}|^{\rm n}\right]^{\rm m}} \tag{3}$$

where Θ is the effective saturation (dimensionless), θ(ψ <sup>m</sup>) is the volumetric soil moisture (m3/m<sup>3</sup> ), θr is the residual volumetric soil moisture (m3/m<sup>3</sup> ), θs is the volumetric soil moisture at saturation (m3/m<sup>3</sup> ), ψ<sup>m</sup> is the matric potential of water in the soil (kPa), α the value of air input per kPa, n is the empirical parameter of the fit (dimensionless), and m the Burdine restriction where m = 1 − (2/n). The volumetric soil moisture (θhco) measured in m3/m<sup>3</sup> at the pore water suction point in the hydraulic cut-off (ψhco measured in kPa) was calculated according to Czyz and Dexter 2013 [25] considering ψhco = −1030 kPa.

#### *2.5. Unsaturated Soil Hydraulic Conductivity and Soil Water Redistribution*

The unsaturated soil hydraulic conductivity (K(θ)), measured in cm/hour, was calculated using the van Genuchten–Burdine equation:

$$\mathbf{K}(\boldsymbol{\Theta}) = \mathbf{K}\_{\mathbf{S}} \cdot \boldsymbol{\Theta}^{2} \cdot \left[ 1 - \left( 1 - \boldsymbol{\Theta}^{\frac{1}{\mathbf{m}}} \right)^{\mathbf{m}} \right] \tag{4}$$

The volumetric soil moisture (θ) at potentials −1 to −6 kPa was calculated (for matric potential or ψ less than −6 kPa the water is retained) to estimate K(θ). Subsequently, the K(ψ) plot was prepared for soil layers at 0 to 0.2 m, 0.2 to 0.4 m, and 0.4 to 0.6 m depths.

In order to estimate the effective saturation (Θ), which was calculated using the actual volumetric moisture at field capacity (θˆ CC):

$$\Theta = \frac{\mathfrak{G}\_{\text{CC}} - \mathfrak{Re}}{\mathfrak{G}\mathfrak{s} - \mathfrak{H}\mathfrak{r}} = \frac{1}{\left(1 + \left|\mathfrak{a} \cdot \psi\_{\text{CC}}\right|^{\mathfrak{n}}\right)^{\mathfrak{m}}} \tag{5}$$

and in order to estimate the matric potential at field capacity (ψCC):

$$\psi\_{\text{CC}} = \frac{\left[\left(\frac{\Theta \text{s} - \Theta \text{r}}{\Theta\_{\text{CC}} - \Theta \text{r}}\right)\frac{1}{\text{m}} - 1\right]^{\frac{1}{\text{n}}}}{\alpha} \tag{6}$$

It was assumed that this condition was reached when the drainage rate (τ) was defined as a percentage (p) of the saturated hydraulic conductivity (KS) based on Prevedello 1999 [26] as restricted by Burdine 1953 [23]. The equations involving this assumption are as follows:

$$
\boldsymbol{\pi} = \mathbf{K}(\boldsymbol{\theta}) \tag{7}
$$

$$\boldsymbol{\pi} = \mathbf{K}\_{\rm S} \cdot \boldsymbol{\Theta}^{\rm 2} \cdot \left[ 1 - \left( 1 - \boldsymbol{\Theta}^{\frac{1}{\rm m}} \right)^{\rm m} \right] \tag{8}$$

$$\mathbf{p} = \frac{\mathbf{r}}{\mathbf{K}\_{\mathbf{S}}} = \boldsymbol{\Theta}^{\frac{\mathbf{r}}{2}} \cdot \left[ 1 - \left( 1 - \boldsymbol{\Theta}^{\frac{1}{\mathbf{m}}} \right)^{\mathbf{m}} \right] \tag{9}$$

Based on Sisson et al., 1980 (Equation (10)) [27] and Prevedello and Armindo 2015 (Equation (11)) [28], an analytical model of soil water redistribution was obtained. For this, Equation (4) was derived with respect to θ and equalized to z/t:

$$\frac{\mathbf{z}}{\mathbf{t}} = \frac{\mathbf{d}\mathbf{K}}{\mathbf{d}\theta} \tag{10}$$

$$\frac{d\mathbf{K}}{d\theta} = \frac{d\mathbf{K}}{d\Theta} \cdot \frac{d\Theta}{d\theta} \tag{11}$$

Then, the time corresponding to the field capacity in each soil layer was obtained:

$$\frac{\mathbf{z}}{\mathbf{f}\_{\mathbf{C}\mathbf{C}}} = \left[ \left( \mathbf{K}\_{\mathbf{S}} \cdot (\mathbf{2} \cdot \boldsymbol{\Theta}) \cdot \left( 1 - \left( 1 - (\boldsymbol{\Theta})^{\frac{1}{\mathbf{m}}} \right)^{\mathbf{m}} \right) + \mathbf{K}\_{\mathbf{S}} \cdot (\boldsymbol{\Theta})^2 \cdot \left( \left( 1 - (\boldsymbol{\Theta})^{\frac{1}{\mathbf{m}}} \right)^{\mathbf{m}-1} \cdot \left( \mathbf{m} \cdot \left( (\boldsymbol{\Theta})^{\frac{1}{\mathbf{m}} - 1} \cdot \frac{1}{\mathbf{m}} \right) \right) \right) \right) \cdot \left( \frac{1}{(\boldsymbol{\Theta} \mathbf{s} - \boldsymbol{\theta} \mathbf{r})} \right) \right] \tag{12}$$

where p is the percentage of K<sup>S</sup> (where we considered the values: 0.005; 0.010; 0.015; 0.020; 0.025; 0.030; 0.040; and 0.050), Θ is the effective saturation (dimensionless), m is the parameter from the Burdine (1953) restriction (dimensionless) [23]; tCC is the time corresponding to the field capacity (h), z is the soil layer (m), and K<sup>S</sup> is the saturated hydraulic conductivity (m/h). Estimation of the Θ function was performed in a routine developed in the R statistical program [24]. The moisture at the inflection point (θˆ IP measured in m3/m<sup>3</sup> ) of the soil water retention curve was determined according to Dexter and Bird (2001) [29], and the respective matric potential at the inflection point (ψIP measured in kPa) was estimated using Equation (9) as:

$$\hat{\boldsymbol{\Theta}}\_{\rm IP} = \begin{pmatrix} \boldsymbol{\Theta}\mathbf{s} \ - & \boldsymbol{\Theta}\mathbf{r} \end{pmatrix} \cdot \left(\mathbf{1} + \frac{1}{\mathbf{m}}\right)^{-\mathbf{m}} + \boldsymbol{\Theta}\mathbf{r} \tag{13}$$

#### *2.6. Soil Water Capacity Function*

The soil pore size distribution, represented by the water capacity function, was estimated using the retention curve model and capillarity theory [3,30]. The van Genuchten (1980) model was modified by replacing the matric potential ψ<sup>m</sup> with the pore radius (r), using the capillarity equation [3]:

$$\boldsymbol{\Theta} = \boldsymbol{\theta}\mathbf{r} + \frac{\left(\boldsymbol{\theta}\mathbf{s} - \,\,\theta\mathbf{r}\right)}{\left[1 + \left(\frac{\mathbf{A}}{\mathbf{r}}\right)^{\mathbf{n}}\right]^{\mathbf{m}}}\tag{14}$$

$$\mathbf{A} = \mathbf{2} \cdot \boldsymbol{\sigma} \cdot \boldsymbol{\alpha} \cdot 10^3 \tag{15}$$

With the derivative of the modified model, one obtains the following expressions for rmax (maximum pore radius) and dmax (maximum pore diameter):

$$\frac{\text{d}\theta}{\text{d}\text{lorg}} = (\text{\(\theta s - \theta r\)} \cdot \text{m} \cdot \text{n} \cdot \text{A}^{\text{n}} \cdot \text{r}^{-\text{n}} \cdot (1 + \text{A}^{\text{n}} \cdot \text{r}^{-\text{n}})^{-\text{m}-1} \tag{16}$$

Deriving and equating to zero, one has:

$$\mathbf{r}\_{\text{max}} = \mathbf{A} \cdot \left(\frac{1}{\mathbf{m}}\right)^{-\frac{1}{\mathbf{n}}} \tag{17}$$

$$\mathbf{d}\_{\text{max}} = \mathbf{2} \cdot \mathbf{r}\_{\text{max}} \tag{18}$$

where σ is the surface tension of water (N m−<sup>1</sup> ); r is the pore radius measured in micrometers (µm).

#### *2.7. Theoretical and Actual Soil Water Storage and Availability*

The available soil water capacity (AWC) was determined from soil moisture measured at −6 kPa (AWC1) and estimated from matric potential at field capacity (AWC2), relative to the −1500 kPa matric potential:

$$\text{AWC}\_{1} = \sum\_{i=1}^{n} (\theta\_{\text{CCi}} - \theta\_{\text{PMPi}}) \cdot \text{z}\_{i} \tag{19}$$

$$\text{AWC}\_2 = \sum\_{i=1}^{n} \left(\hat{\theta}\_{\text{CCi}} - \theta\_{\text{PMPi}}\right) \cdot \text{z}\_i \tag{20}$$

where AWC is measured in millimeters; θCCi the volumetric soil moisture at field capacity in the *i*th layer (m3/m<sup>3</sup> ) at the matric potential of −6 kPa; θPMPi the volumetric soil moisture at the permanent wilting point in the *i*th layer (m3/m<sup>3</sup> ) at the matric potential of −1500 kPa, θˆ CC the estimated soil volumetric moisture at field capacity in the *i*th layer (m3/m<sup>3</sup> ) at the matric potential ψCC; z<sup>i</sup> is the thickness of the *i*th soil layer (mm); and n the number of layers considered, which is dimensionless.

Considering the results obtained from the soil water redistribution and point of hydraulic cut-off, the easily available water (EAW) at matric potential −6 kPa (EAW1) and ψCC (EAW2) was determined in relation to the volumetric moisture at 1030 kPa −1030 kPa (ψhco):

$$\text{EAW}\_1 = \sum\_{\mathbf{i}=1}^{n} (\theta\_{\text{CCi}} - \theta\_{\text{hco}}) \cdot \mathbf{z}\_{\mathbf{i}} \tag{21}$$

$$\text{EAW}\_2 = \sum\_{i=1}^{n} \left(\hat{\theta}\_{\text{CCi}} - \theta\_{\text{hco}}\right) \cdot \mathbf{z}\_{\text{i}} \tag{22}$$

where EAW in the soil is measured in mm, θhco is the volumetric soil moisture at the matric potential of <sup>−</sup>1030 kPa in the *<sup>i</sup>*th layer (m3/m<sup>3</sup> ), z<sup>i</sup> is the thickness of the *i*th soil layer (mm), and n is the number of layers considered and is dimensionless. Available soil water capacity (AWC) and easily available water (EAW) were calculated for different soil layers at four percentages of AWC and EAW: 100%, 80%, 60%, and 40%. These values were then summarized using R software's *ggplot2* function [24].

#### *2.8. Statistical Analyses*

Descriptive statistics were used to analyze our results with the help of the boxplot graph, and the discriminant values (outliers) were removed from the subsequent analyses with 6 repetitions. The assumptions of normality and homogeneity of variances were verified with the Shapiro–Wilk and Bartlett tests, respectively. Once these assumptions were met, the data were submitted for analysis of variance and the averages were compared using the Tukey test at the 5% significance level. The tCC data were submitted to the Kruskal– Wallis test. Statistical errors and index between observed (Y) and estimated (Yˆ ) values were quantified with the following expressions for root mean square error (RMSE), the root mean square error normalized by the mean (NRMSE), the ratio of RMSE to the standard deviation (RSR), and Willmott's agreement index (d) [31,32]:

$$\text{RMSE} = \sqrt{\frac{\sum\_{i=1}^{n} \left(\hat{\mathbf{Y}}\_{i} - \mathbf{Y}\_{i}\right)^{2}}{n}} \tag{23}$$

$$\text{NRMSE} = \frac{\sqrt{\frac{\sum\_{i=1}^{n} \left(\mathbf{\hat{Y}}\_{i} - \mathbf{Y}\_{i}\right)^{2}}{n}}}{\overline{\mathbf{Y}}} \tag{24}$$

$$\text{RSR} = \frac{\text{RMSE}}{\text{Dp}\_{\text{Y}}} \tag{25}$$

$$\mathbf{d} = 1 - \left[ \frac{\sum\_{i=1}^{n} \left( \mathbf{\hat{Y}\_i} - \mathbf{Y\_i} \right)^2}{\sum\_{i=1}^{n} \left( \left| \mathbf{\hat{Y}\_i} - \overline{\mathbf{Y}} \right| + \left| \mathbf{Y\_i} - \overline{\mathbf{Y}} \right| \right)^2} \right] \tag{26}$$

with measurement units specified for RMSE (m3/m<sup>3</sup> ), NRMSE (%), RSR (dimensionless), and d (dimensionless). For RMSE, NRMSE, and d, Yˆ i is the *i*th value of the estimated variable (m3/m<sup>3</sup> ), Y<sup>i</sup> is the *i*th value of the observed variable (m3/m<sup>3</sup> ), and n is the amount of data for the observed variable (units of variable). When calculating NRMSE and d, Y is the mean of the values of the observed variable (m3/m<sup>3</sup> ). For RSR calculation, Dp<sup>Y</sup> is the standard deviation of the observed data (m3/m<sup>3</sup> ).

Validation of the associations among soil physical-hydric attributes was performed using principal component analysis (PCA). The correlation matrix of the standardized variables (unit variances and zero means) of each soil layer was used in PCA. The *prcomp* function was used to perform PCA and the *ggfortify* package [33] of the R software [24] was used for graphing.

#### *2.9. Sugarcane Yield Modeling*

The 12 years at the research site where sugarcane cultivars were evaluated were 1998–2006, 2008, and 2018–2019. Over these years, there were four sugarcane cultivars that were tested. RB72454 was planted from 1998 to 2006 and in 2008. RB867515 was trialed from 2004 to 2006 and from 2018 to 2019. RB966928 was evaluated from 2004 to 2005 and from 2018 to 2019. Finally, RB036066 was planted in 2018 and 2019. Sugarcane was grown in all years without irrigation.

Simulated sugarcane yields were estimated using soil water storage values (AWC and EAW, considered at 100%) from this study entered into a model validated from 12 years of field data from UFPR and RIDESA by Viana et al., 2023 [34]:

$$\text{TCH} = -986.016 - 0.288 \times \text{ADD}\_{\text{I}} + 165.645 \times \text{SWS}\_{\text{II}} \tag{27}$$

In the simulation, the logarithmic transformation of the multiple linear regression model was performed:

$$\text{TCH}=\text{e}^{-11.3374} \times \text{ADD}\_{\text{I}}{}^{-0.8224} \times \text{SWS}\_{\text{II}}{}^{10.5272} \tag{28}$$

where TCH is metric tons of sugarcane stalk per hectare, ADD<sup>I</sup> is growing degree days ( ◦C) during the first stage of sugarcane crop development (July through October), and SWSII is soil water storage during the second development (October through March) for a 12-month production period. The third stage of development runs from March through July where the crop accumulates sucrose [34]. Simulated yields (*y*-axis) were then graphed three-dimensionally against SWSII (AWC and EAW) measured in cm (*x*-axis) and ADD<sup>I</sup> measured in ◦C (*z*-axis) for all data points. Graphing was performed using R's *scatter3d plotly r* function [24].

#### **3. Results**

#### *3.1. Texture and Sandy Soil Structure*

Clay content was significantly higher in the 0.2 to 0.4 and 0.4 to 0.6 m layers of soil compared to the topsoil layer at 0 to 0.2 m depth (Table 1). There was an inverse relationship between coarse and fine sand fractions in the soil profile (Table 1). The soil bulk density (ρS) and particle density (ρP) did not show significant differences among the soil layers (Table 1). The magnitude of the values we found is consistent with what is expected from sandy texture soils [1]. The soil bulk densities we measured (Table 1) ranged from 1582 to 1690 kg/m<sup>3</sup> which were consistent with an expected range of 1300 to 1800 kg/m<sup>3</sup> [30]. The average value of total porosity (TP = 0.33 m3/m<sup>3</sup> ) was low, as commonly observed in sandy soils with fine sand contents [2,35]. As macroporosity decreased, microporosity increased, but ρ<sup>S</sup> did not increase and TP decreased (*p* < 0.05) in the subsurface layers (Table 1). There was an increase in microporosity and clay contents at depth.


**Table 1.** Physical-hydric attributes of *Latossolo Vermelho distrófico* (Oxisol) soil layers at the Sugarcane Genetic Improvement Program Experimental Station cultivated with sugarcane in Paranavaí, Paraná state, Brazil. **Table 1.** Physical-hydric attributes of *Latossolo Vermelho distrófico* (Oxisol) soil layers at the Sugarcane Genetic Improvement Program Experimental Station cultivated with sugarcane in Paranavaí, Paraná state, Brazil.

<sup>1</sup> Averages followed by equal letters in the column do not differ by the Tukey test at 5% probability. <sup>2</sup> Var is the variance (variable unit). <sup>3</sup> Stdev is the standard deviation (variable unit). <sup>4</sup> CV is the coefficient of variation (%). coefficient of variation (%).

#### *3.2. Water Movement in Sandy Soil 3.2. Water Movement in Sandy Soil*

The saturated hydraulic conductivity (KS) was higher in the 0 to 0.2 m layer and reduced abruptly at depth (Figure 1a) due to macroporosity facilitating water drainage in saturated soil [35]. Like Ks, unsaturated hydraulic conductivity K(θ) was also higher in the uppermost soil layer. However, the decrease in K(θ) with the reduction in soil moisture was more significant in the 0 to 0.2 m layer, compared to the others (Figure 1b). The K(ψ) in the surface layer showed greater decay (Figure 1c) with increasing matric potential due to the smaller number of capillary pores [35]. The saturated hydraulic conductivity (KS) was higher in the 0 to 0.2 m layer and reduced abruptly at depth (Figure 1a) due to macroporosity facilitating water drainage in saturated soil [35]. Like Ks, unsaturated hydraulic conductivity K(θ) was also higher in the uppermost soil layer. However, the decrease in K(θ) with the reduction in soil moisture was more significant in the 0 to 0.2 m layer, compared to the others (Figure 1b). The K(ψ) in the surface layer showed greater decay (Figure 1c) with increasing matric potential due to the smaller number of capillary pores [35].

**Figure 1***.* (**a**) Saturated (KS) conductivity with error bars in cm/hour; (**b**) unsaturated (K(θ)) hydraulic conductivity versus volumetric water content; and (**c**) versus matric potential in the soil layers 0 to 0.20 m (m), 0.2 to 0.4 m, and 0.4 to 0.6 m for the Latossolo Vermelho distrófico (Oxisol) at the Sugarcane Genetic Improvement Program Experimental Station cultivated with sugarcane in Paranavaí, Paraná state, Brazil. **Figure 1.** (**a**) Saturated (KS) conductivity with error bars in cm/hour; (**b**) unsaturated (K(θ)) hydraulic conductivity versus volumetric water content; and (**c**) versus matric potential in the soil layers 0 to 0.20 m (m), 0.2 to 0.4 m, and 0.4 to 0.6 m for the Latossolo Vermelho distrófico (Oxisol) at the Sugarcane Genetic Improvement Program Experimental Station cultivated with sugarcane in Paranavaí, Paraná state, Brazil.

#### *3.3. Water Retention and Pore Distribution 3.3. Water Retention and Pore Distribution*

The soil water retention curves showed no evidence of compaction effects (Figure 2), as the soil bulk density (ρS) values are below the critical value of ρ<sup>S</sup> > 1.70 g/cm3 [2]. The parameters of the van Genuchten model (α and n) were significant (p < 0.05) to represent the soil water retention curve for all soil layers (Table 2). Statistical errors and indices between observed (θ) and estimated (θ) soil volumetric moisture indicated that fits with the van Genuchten–Burdine model were satisfactory for the analyzed layers (Table 3). Among The soil water retention curves showed no evidence of compaction effects (Figure 2), as the soil bulk density (ρS) values are below the critical value of ρ<sup>S</sup> > 1.70 g/cm<sup>3</sup> [2]. The parameters of the van Genuchten model (α and n) were significant (*p* < 0.05) to represent the soil water retention curve for all soil layers (Table 2). Statistical errors and indices between observed (θ) and estimated (θˆ) soil volumetric moisture indicated that fits with the van Genuchten–Burdine model were satisfactory for the analyzed layers (Table 3). Among the errors, the NRMSE was higher due to the use of the six repetitions at each matric potential to perform the adjustment of the soil water retention curve, which improved the significance of the adjustment.

significance of the adjustment.

**Figure 2.** The water retention curve of the Latossolo Vermelho distrófico (Oxisol) in the layers (**a**) 0 to 0.2 m; (**b**) 0.2 to 0.4 m; and (**c**) 0.4 to 0.6 m. The fitted values are accompanied by the confidence band for the predicted value. **Figure 2.** The water retention curve of the Latossolo Vermelho distrófico (Oxisol) in the layers (**a**) 0 to 0.2 m; (**b**) 0.2 to 0.4 m; and (**c**) 0.4 to 0.6 m. The fitted values are accompanied by the confidence band for the predicted value.

The RSR is one of the quantitative statistics recommended by Moriasi et al. (2007) [36] to evaluate models. The smallest magnitudes of RSR were observed at different p in the soil layers due to the physical-hydric attributes of the soil (Table 4). Lower matric potential ψIP was found in the 0 to 0.2 m (m) layer due to the soil texture (Table 5). The pore distribution curve (Figure 3) showed a maximum pore diameter (dmax) of 64.19 µm (µm) at a soil depth

of 0 to 0.2 m, 51.96 µm at 0.2 to 0.4 m, and 50.6 µm at 0.4 to 0.6 m. The maximum point of the pore diameter frequency curve dθ/dlog(r) corresponded to 0.192 at a depth of 0 to 0.2 m, 0.1458 at 0.2 to 0.4 m, and 0.1463 at 0.4 to 0.6 m in the analyzed soil layers. Matric potential at field capacity (ψCC in kPa) <sup>1</sup> 0 to 0.2 6.55 c 5.93 c 5.59 c 5.35 c 5.16 b 5.01 b 4.78 b 4.60 b 0.2 to 0.4 8.16 ab 7.38 ab 6.95 ab 6.64 ab 6.41 ab 6.23 ab 5.94 ab 5.71 ab


0.4 to 0.6 8.56 a 7.73 a 7.26 a 6.94 a 6.69 a 6.49 a 6.18 a 5.94 a

**Table 2.** Fitted parameters for the van Genuchten model and their confidence intervals. Average 7.76 7.01 6.60 6.31 6.09 5.91 5.63 5.42

> <sup>1</sup> θs and θr are the soil moisture at saturation and residual, respectively, while α and n are the parameters of the van Genuchten model. **Soil Layer Estimated Moisture Matric Potential <sup>ψ</sup>IP (kPa) <sup>1</sup> RSR (Dimensionless) <sup>5</sup>**

> **Table 3.** Errors and statistical index between observed (θ) and estimated (θˆ) soil volumetric moisture with van Genuchten model using Burdine restriction. 0 to 0.2 0.223 b 4.48 b 4.29 0.2 to 0.4 0.234 ab 5.37 ab 1.08 0.4 to 0.6 0.244 a 5.79 a 1.37


**(m3/m3) <sup>1</sup>**

**Depth (meter)**

**Figure 3.** Pore distribution curve for soil layers at depths of 0 to 0.2 m, 0.2 to 0.4 m*,* and 0.4 to 0.6 m of Latossolo Vermelho distrófico (Oxisol) used for sugarcane cultivation at the Experimental Station of the Sugarcane Genetic Improvement Program at the Federal University of Paraná and the Inter-University Network for the Development of Sugarcane Energy. **Figure 3.** Pore distribution curve for soil layers at depths of 0 to 0.2 m, 0.2 to 0.4 m, and 0.4 to 0.6 m of Latossolo Vermelho distrófico (Oxisol) used for sugarcane cultivation at the Experimental Station of the Sugarcane Genetic Improvement Program at the Federal University of Paraná and the Inter-University Network for the Development of Sugarcane Energy.


**Table 4.** Estimated actual moisture (θˆ CC), matric potential (ψCC), and time (tCC) at field capacity and statistical error between observed (θCC) at –6 kPa and estimated (θˆ CC) volumetric moisture at field capacity.

<sup>1</sup> Averages followed by equal letters in the column do not differ by Tukey's test at 5% probability (θˆ CC and ψCC) and tCC by Kruskal–Wallis test.


**Table 5.** Estimated moisture (θˆ IP) and matric potential (ψIP) at the inflection point, and statistical error between observed (θCC) at –6 kPa and estimated (θˆ IP) volumetric moisture at field capacity.

<sup>1</sup> Averages followed by equal letters in the column do not differ by the Tukey test at 5% probability. <sup>2</sup> Var is the variance (variable unit). <sup>3</sup> Stdev is the standard deviation (variable unit). <sup>4</sup> CV is the coefficient of variation (%). <sup>5</sup> RSR—Root mean square error/Standard deviation (dimensionless).

#### *3.4. Water Availability*

In general, the soil of the experimental area cultivated with sugarcane showed low water availability for plants in the layers analyzed (Table 6), which is characteristic of soils with low clay contents [5,37]. There was no statistical difference between the layers (*p* > 0.05), but the highest AWC<sup>1</sup> and EAW<sup>1</sup> occur in the soil layers at 0.2 to 0.4 m and 0.4 to 0.6 m depth due to the higher amount of clay and microporosity [5]. The highest AWC<sup>2</sup> and EAW<sup>2</sup> occurred in the surface layer (*p* < 0.05) due to the higher organic carbon content that increases water retention and availability for plants [38]. This surface layer corresponds to the soil layer with the highest volume and density of root length for sugarcane [34]. Figure 4 shows the variation in the percentage of water availability for the soil in the experimental area at the Experimental Station of the Sugarcane Genetic Improvement Program at the Federal University of Paraná and the Inter-University Network for the Development of Sugarcane Energy. 0.05), but the highest AWC1 and EAW1 occur in the soil layers at 0.2 to 0.4 m and 0.4 to 0.6 m depth due to the higher amount of clay and microporosity [5]. The highest AWC2 and EAW2 occurred in the surface layer (*p* < 0.05) due to the higher organic carbon content that increases water retention and availability for plants [38]. This surface layer corresponds to the soil layer with the highest volume and density of root length for sugarcane [34]. Figure 4 shows the variation in the percentage of water availability for the soil in the experimental area at the Experimental Station of the Sugarcane Genetic Improvement Program at the Federal University of Paraná and the Inter-University Network for the Development of Sugarcane Energy.

In general, the soil of the experimental area cultivated with sugarcane showed low water availability for plants in the layers analyzed (Table 6), which is characteristic of soils with low clay contents [5,37]. There was no statistical difference between the layers (*p* >

**Table 6.** Hydric attributes of the *Latossolo Vermelho distrófico* (Oxisol), cultivated with sugarcane at the Experimental Station of the Sugarcane Genetic Improvement Program at the Federal University of Paraná and the Inter-University Network for the Development of Sugarcane Energy. **Table 6.** Hydric attributes of the *Latossolo Vermelho distrófico* (Oxisol), cultivated with sugarcane at the Experimental Station of the Sugarcane Genetic Improvement Program at the Federal University of Paraná and the Inter-University Network for the Development of Sugarcane Energy.


*Sustainability* **2023**, *15*, x FOR PEER REVIEW 12 of 23

*3.4. Water Availability*

<sup>1</sup> Averages followed by equal letters in the column do not differ by the Tukey test at 5% probability. <sup>2</sup> Var is the variance (variable unit). <sup>3</sup> Stdev is the standard deviation (variable unit). <sup>4</sup> CV is the coefficient of variation (%). <sup>2</sup> Var is the variance (variable unit). <sup>3</sup> Stdev is the standard deviation (variable unit). <sup>4</sup> CV is the coefficient of variation (%).

#### *3.5. Multivariate Relationships among Soil Physical‐Hydric Attributes 3.5. Multivariate Relationships among Soil Physical-Hydric Attributes*

Principal component analysis (PCA) was used to validate the previously established associations between the structural and hydric physical attributes. We observed that the first principal component (PC1) explained 60.6% of the total variance. PC1 (eigenvalue of 3.114) and PC2 (eigenvalue of 1.495) together explained 74.6% of the total variance (Figure 5). Principal component analysis (PCA) was used to validate the previously established associations between the structural and hydric physical attributes. We observed that the first principal component (PC1) explained 60.6% of the total variance. PC1 (eigenvalue of 3.114) and PC2 (eigenvalue of 1.495) together explained 74.6% of the total variance (Figure 5).

PCA was efficient in distinguishing between soil layers. The 0 to 0.2 m layer was characterized and influenced mainly by the physical-hydric attributes of macroporosity, silt, saturated hydraulic conductivity, unsaturated hydraulic conductivity, sand, coarse sand, and maximum pore diameter. The variables fine sand, clay, microporosity, volumetric soil moisture at the permanent wilting point, volumetric soil moisture at the inflection point, soil density, and available soil water capacity determined at –6 kPa were higher and more expressive in the physical-hydric attributes of the soil layers at depths of 0.2 to 0.4 m and 0.4 to 0.6 m (Figure 5). It was found that saturated hydraulic conductivity showed a direct relationship with sand, coarse sand, silt, macroporosity, and maximum pore diameter, and inverse with microporosity, clay, volumetric soil moisture at the permanent wilting point, and fine sand. The soil density showed an inverse relationship with macroporosity.

Microporosity was directly related to clay, fine sand, and volumetric soil moisture at the permanent wilting point, and inversely with saturated hydraulic conductivity, sand, coarse sand, silt, macroporosity, and maximum pore diameter. *Sustainability* **2023**, *15*, x FOR PEER REVIEW 13 of 23

> **Figure 5.** Principal component analysis (PCA) of the physical-hydric attributes of the Latossolo Vermelho distrófico (Oxisol), cultivated with sugarcane in the Experimental Station of the Sugarcane Genetic Improvement Program at the Federal University of Paraná and the Inter-University Network for the Development of Sugarcane Energy. Labels for (1) "dmax" and "Coarse" and (2) "Clay" and "Micro" slightly overlap with each other due to close proximity of these two pairs of vectors. **Figure 5.** Principal component analysis (PCA) of the physical-hydric attributes of the Latossolo Vermelho distrófico (Oxisol), cultivated with sugarcane in the Experimental Station of the Sugarcane Genetic Improvement Program at the Federal University of Paraná and the Inter-University Network for the Development of Sugarcane Energy. Labels for (1) "dmax" and "Coarse" and (2) "Clay" and "Micro" slightly overlap with each other due to close proximity of these two pairs of vectors.

#### PCA was efficient in distinguishing between soil layers. The 0 to 0.2 m layer was *3.6. Sugarcane Yield Simulations*

characterized and influenced mainly by the physical-hydric attributes of macroporosity, silt, saturated hydraulic conductivity, unsaturated hydraulic conductivity, sand, coarse sand, and maximum pore diameter. The variables fine sand, clay, microporosity, volumetric soil moisture at the permanent wilting point, volumetric soil moisture at the inflection point, soil density, and available soil water capacity determined at –6 kPa were higher and more expressive in the physical-hydric attributes of the soil layers at depths of 0.2 to 0.4 m and 0.4 to 0.6 m (Figure 5). It was found that saturated hydraulic conductivity showed a direct relationship with sand, coarse sand, silt, macroporosity, and maximum pore diameter, and inverse with microporosity, clay, volumetric soil moisture at the permanent For 12 years at the research site, sugarcane stalk yield (TCH) averaged 114.68 metric tons (t)/hectare (ha). TCH was at a minimum of 65.997 t/ha in 2005 and a maximum of 164.190 t/ha in 2006 (Figure 6). Fracaro et al., 2014 [39] found sugarcane TCH of 170 t/ha for RB036152 and an average TCH of 110 t/ha across 26 genotypes from RIDESA in the 2013–2014 crop season in the municipality of Viamão, Rio Grande do Sul state, Brazil. The estimated Brazilian average national productivity of sugarcane for the 2022–2023 crop season, according to CONAB [40], is 72 t/ha. The 2022–2023 crop season was characterized by reduced rainfall and low temperatures in the Central-South Region, which accounts for around 90% of Brazil's total sugarcane production [40].

wilting point, and fine sand. The soil density showed an inverse relationship with macroporosity. Microporosity was directly related to clay, fine sand, and volumetric soil moisture at the permanent wilting point, and inversely with saturated hydraulic conductivity, sand, coarse sand, silt, macroporosity, and maximum pore diameter. *3.6. Sugarcane Yield Simulations* For 12 years at the research site, sugarcane stalk yield (TCH) averaged 114.68 metric tons (t)/hectare (ha). TCH was at a minimum of 65.997 t/ha in 2005 and a maximum of 164.190 t/ha in 2006 (Figure 6). Fracaro et al., 2014 [39] found sugarcane TCH of 170 t/ha for RB036152 and an average TCH of 110 t/ha across 26 genotypes from RIDESA in the In the sugarcane model simulation with AWC and EAW set at a value of 100% (Figure 7), TCH for: (1) EAW<sup>1</sup> averaged 33.94 t/ha with a minimum of 26.84 t/ha and a maximum of 42.27 t/ha; (2) AWC<sup>1</sup> averaged 72.10 t/ha with a minimum of 57.00 t/ha and a maximum of 89.78 t/ha; (3) EAW<sup>2</sup> averaged 153.96 t/ha with a minimum of 121.73 t/ha and a maximum of 191.73 t/ha; and (4) AWC<sup>2</sup> averaged 296.58 t/ha with a minimum of 234.50 t/ha and a maximum of 369.34 t/ha. Under the 100% water availability condition in Phase II, EAW<sup>1</sup> underestimates and AWC<sup>2</sup> overestimates sugarcane yield. Analyzing the simulated values with EAW<sup>2</sup> (set at 100%), the potential yield of sugarcane can be achieved with sustainable management practices, genetic improvement, and favorable climate conditions.

2013–2014 crop season in the municipality of Viamão, Rio Grande do Sul state, Brazil. The estimated Brazilian average national productivity of sugarcane for the 2022–2023 crop season, according to CONAB [40], is 72 t/ha. The 2022–2023 crop season was characterized by reduced rainfall and low temperatures in the Central-South Region, which accounts tions.

**Figure 6.** Sugarcane yield model validated from 12 years of field data from UFPR and RIDESA by Viana et al., 2023 [34]. **Figure 6.** Sugarcane yield model validated from 12 years of field data from UFPR and RIDESA by Viana et al., 2023 [34].

In the sugarcane model simulation with AWC and EAW set at a value of 100% (Figure 7), TCH for: (1) EAW1 averaged 33.94 t/ha with a minimum of 26.84 t/ha and a maximum of 42.27 t/ha; (2) AWC1 averaged 72.10 t/ha with a minimum of 57.00 t/ha and a maximum of 89.78 t/ha; (3) EAW2 averaged 153.96 t/ha with a minimum of 121.73 t/ha and a maximum of 191.73 t/ha; and (4) AWC2 averaged 296.58 t/ha with a minimum of 234.50 t/ha and a maximum of 369.34 t/ha. Under the 100% water availability condition in Phase II, EAW1 underestimates and AWC2 overestimates sugarcane yield. Analyzing the simulated values with EAW2 (set at 100%), the potential yield of sugarcane can be achieved with sustainable management practices, genetic improvement, and favorable climate condi-

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 15 of 23

**Figure 7.** Simulated sugarcane stalk yield (TCH = t/ha) as a function of growing degree days during development phase I (ADDI) and soil water storage during development phase II (SWSII) using both AWC and EAW (100%), with the logarithmic transformation of the validated model by Viana et al., 2023 [34]. **Figure 7.** Simulated sugarcane stalk yield (TCH = t/ha) as a function of growing degree days during development phase I (ADD<sup>I</sup> ) and soil water storage during development phase II (SWSII) using both AWC and EAW (100%), with the logarithmic transformation of the validated model by Viana et al., 2023 [34].

#### **4. Discussion 4. Discussion**

*4.1. Comparisons to Previous Research 4.1. Comparisons to Previous Research*

adopted to minimize soil compaction.

#### 4.1.1. Soil Characteristics 4.1.1. Soil Characteristics

Although there was an increase in clay content in the subsurface layers, this did not result in abrupt textural changes in the soil. Greater clay content in lower soil surface lay-Although there was an increase in clay content in the subsurface layers, this did not result in abrupt textural changes in the soil. Greater clay content in lower soil surface layers

ers can result in increased water retention for sandy soils. As a consequence, the sand and silt contents were lower in the 0.2 to 0.4 m and 0.4 to 0.6 m soil layers. Thus, the analyzed

Sandstone Formation in the northwestern region of Paraná [6]. The higher coarse sand content of the sand fraction is characteristic of soils of the Paranavaí region, confirmed by the granulometric characterization in areas of crop-livestock and pineapple integration in the Caiuá and Paranavaí Sandstone Formations [6]. Fidalski 2017 [42] considered that the granulometric characterization is a criterion to guide the use and management of soils of the Caiuá Sandstone in the northwestern region of Paraná. This is due to the inverse rela-

Souza et al., 2015 [37] found an average soil bulk density of 1620 kg/m3 at the RIDESA. This was close to the average ρ<sup>S</sup> obtained in this study of 1640 kg/m3. In the literature, there are reports of root growth restriction in sugarcane for ρ<sup>S</sup> = 1700 kg/m3 [2]. Thus, the bulk density values obtained in this study between 1600 and 1700 kg/m3 are at the threshold of possibly limiting crop root growth. Therefore, management practices should be

Ortiz et al., 2017 [43] found TP between 0.29 and 0.38 m3/m3 in sandy soil grown in the first, third, and fifth sugarcane cycles in 44, 40, and 15 years of production, respectively. Low macroporosity pore values (0.09 m3/m3) were obtained in the soil layers at 0.2 to 0.4 m and 0.4 to 0.6 m (Table 1). This may be a problem since values lower than 0.1 m3/m3 are limiting to soil aeration and harmful to agricultural production [1,44,45]. These

tionship of the sand contents (total or coarse) with the available water in the soil.

can result in increased water retention for sandy soils. As a consequence, the sand and silt contents were lower in the 0.2 to 0.4 m and 0.4 to 0.6 m soil layers. Thus, the analyzed soil is in the textural classes of sand in the soil layer at the depth range of 0 to 0.20 m, and loamy sand in soil layers 0.2 to 0.4 m and 0.4 to 0.6 m [41], characteristic of the Paranavaí Sandstone Formation in the northwestern region of Paraná [6]. The higher coarse sand content of the sand fraction is characteristic of soils of the Paranavaí region, confirmed by the granulometric characterization in areas of crop-livestock and pineapple integration in the Caiuá and Paranavaí Sandstone Formations [6]. Fidalski 2017 [42] considered that the granulometric characterization is a criterion to guide the use and management of soils of the Caiuá Sandstone in the northwestern region of Paraná. This is due to the inverse relationship of the sand contents (total or coarse) with the available water in the soil.

Souza et al., 2015 [37] found an average soil bulk density of 1620 kg/m<sup>3</sup> at the RIDESA. This was close to the average ρ<sup>S</sup> obtained in this study of 1640 kg/m<sup>3</sup> . In the literature, there are reports of root growth restriction in sugarcane for ρ<sup>S</sup> = 1700 kg/m<sup>3</sup> [2]. Thus, the bulk density values obtained in this study between 1600 and 1700 kg/m<sup>3</sup> are at the threshold of possibly limiting crop root growth. Therefore, management practices should be adopted to minimize soil compaction.

Ortiz et al., 2017 [43] found TP between 0.29 and 0.38 m3/m<sup>3</sup> in sandy soil grown in the first, third, and fifth sugarcane cycles in 44, 40, and 15 years of production, respectively. Low macroporosity pore values (0.09 m3/m<sup>3</sup> ) were obtained in the soil layers at 0.2 to 0.4 m and 0.4 to 0.6 m (Table 1). This may be a problem since values lower than 0.1 m3/m<sup>3</sup> are limiting to soil aeration and harmful to agricultural production [1,44,45]. These results support recommendations for better managing sandy soils, including reducing soil bulk density and increasing total porosity and macroporosity, especially in deeper soil layers [44].

#### 4.1.2. Saturated Soil Hydraulic Conductivity

Our results indicate that the water dynamics in the soils that were analyzed can be modified between their layers since the microporosity is responsible for the retention and conduction of water in unsaturated soil, and macroporosity is responsible for the drainage and aeration of the soil [1,35]. High K<sup>S</sup> values in more superficial layers are characteristic of sandy textured and well-drained soils [5]. However, especially in the 0.4 to 0.6 m layer, the abrupt change in K<sup>S</sup> can make the soil susceptible to shear which may lead to erosion. The combination of high K<sup>S</sup> and sand content account for more than 80% of the variability in soil loss [46]. Our results reinforce the need to increase macroporosity at depth to better manage sugarcane.

Even though the K<sup>S</sup> in the surface layer was higher than the others (Figure 1a), the values were lower than in soils of other textures also cultivated with sugarcane [47,48]. Due to the effects of soil preparation recently performed for sugarcane re-planting such as chiseling and disking, Cherubin et al., 2016 [47] found K<sup>S</sup> (measured in the field) of 250 cm/h in soil with medium clay texture, in the surface layer. The authors found that K<sup>S</sup> showed a direct correlation with total porosity and macroporosity, and an inverse relationship with ρ<sup>S</sup> and the degree of soil compaction, both of which may be directly influenced by soil management. Silva et al., 2018 [48] observed in sandy soils that K<sup>S</sup> decreased with increasing depth in sandy soils, associated with increasing bulk density and microporosity. Silva et al., 2018 [48] found a direct correlation between K<sup>S</sup> with coarse sand and medium sand fractions and an inverse correlation with the available water capacity. In addition, it was found that the available water capacity showed an inverse linear correlation with the coarse sand and medium sand grain size fractions, and direct linear correlation with the micropore volume.

The K(ψ) between soil layers was similar at matric potentials greater than −5 kPa, suggesting that the field capacity of the sandy soil was reached before −6 kPa [26]. The dependent relationship of unsaturated hydraulic conductivity K(θ) with the matric potential indicated that there are modifications in the soil water retention as a consequence of some variation in the pore distribution curve. Thus, the results of K(ψ) of the analyzed sandy soils

demonstrated that the evaluation of soil quality cannot be limited to the separation of TP into macroporosity and microporosity, but a detailed determination of the pore distribution curve is necessary [9].

#### 4.1.3. Soil Water Retention

Higher soil water retention was observed in the subsurface layers (Table 2; Figure 2) due to the increase in microporosity, clay and silt contents, and reduction in coarse sand content at depth (Table 1). The relative error corresponded to 33.5% and 37.6% in the 0.2 to 0.4 m and 0.4 to 0.6 m layers relative to the subsurface layer, respectively. At potentials of −33 kPa, −100 kPa, and −500 kPa, the average relative errors were 35.6% and 39.5%. The largest relative errors (40.3% and 43.9%) occurred at the −1500 kPa matric potential.

Soil retains water due to capillary and adsorption matrix forces, which are responsible for water retention in the capillary pores of aggregates and on the surfaces of soil particles, respectively [49]. The surfaces of clays are negatively charged and adsorb a large number of water molecules on the soil surface due to their high specific surface area. For example, kaolinite can absorb 70,000 to 300,000 cm<sup>2</sup> of water/gram (g) per gram of kaolinite [3,50]. Thus, soil layers with higher clay contents retain more water in the soil, reinforcing the results obtained in this study.

The soils of the Caiuá Sandstone Formation have low clay content with a predominance of kaolinite (>80%), the presence of illite, a ratio of 2:1 clay to minerals, and small amounts of iron and aluminum oxides [16]. The Caiuá Sandstone Formation soils also have a high sand content, with a predominance of large pore diameters that are greater than 30 µm (Figure 3), which is drained at low matric potentials [9,51], reducing the water holding capacity of the soil.

The soil layer with higher fine sand content retained more soil water, especially at low matric potentials (less than −33 kPa; Figure 2). The obtained result reinforces the importance of the fractionation of total sand in sandy soils due to the higher specific surface area of fine sand at 318 cm2/g compared to coarse sand at 79 cm2/g [6,50]. The distribution between fine and coarse sand fractions is an inherent characteristic of soil formation and is not altered by management. Thus, it is worth emphasizing the importance that the accumulation of organic matter has in sandy soils because, besides increasing aggregation, it improves water retention and availability in soils [1,16,38,52,53].

Among the results obtained in the analyses of the soil water retention curves, it was verified that (1) the residual moisture (θr) was similar to the soil volumetric moisture at the hydraulic cut-off point (θhco) in the soil layers at depth ranges of 0 to 0.2 m (0.074 m3/m<sup>3</sup> ), 0.2 to 0.4 m (0.120 m3/m<sup>3</sup> ), and 0.4 to 0.6 m (0.126 m3/m<sup>3</sup> ) (Table 2). Moreover, (2) the residual moisture (θr) was greater than the observed volumetric moisture at –1500 kPa. Finally, (3) the average relative difference in volumetric moisture at –500 and –1500 kPa was 12.92%, 6.03%, and 5.09% in the 0 to 0.2 m, 0.2 to 0.4 m, and 0.4 to 0.6 m layers, respectively. There are indications that the matric potential at the permanent wilting point occurs between −500 kPa and −1030 kPa at the hydraulic cut-off point.

Obtaining the true value of matric potential for improved agricultural crops has been a recurring discussion [54,55]. Torres et al., 2021 [55] found values for matric potential ψPMP < −1500 kPa and found that sunflower, corn, and soybean crops wilted at potentials lower than the estimated matric potentials in the hydraulic cut-off. When plants wilt at water suction lower than the hydraulic cut-off, the permanent wilting point is plant-limited, and when plants wilt at suctions greater than the hydraulic cut-off the limitation is from the soil [25].

Therefore, in sandy soils, obtaining the actual permanent wilting point should receive more attention. In the literature, in sandy texture soil (803 g/kg), there are indications that the average matric potential at the physiological permanent wilting point was –163.7 kPa (Mestre) and –241.7 kPa (Noble) for two modern wheat cultivars, and –90.9 kPa (ANAG 01) and –45.9 kPa (BRS Brau) for two barley cultivars [54].

The relative difference in volumetric moisture at –200 and –1030 kPa was 0.164%, 0.125% and 0.166% in soil layers 0 to 0.2 m, 0.2 to 0.4 m, and 0.4 to 0.6 m, respectively. However, it is critical to emphasize that no studies were found in the literature reporting the estimated matric potential at the physiological permanent wilting point for modern sugarcane cultivars from breeding programs. The hydraulic cut-off suction (–1030 kPa) lower than –1500 kPa reinforces the need for further studies for the sandy soil type analyzed.

Prevedello 1999 [26] found, in homogeneous and heterogeneous soil, that the drainage rate (τ) corresponded to *p* = 0.01 of the saturated soil hydraulic conductivity (τ = 0.01, KS). Andrade and Stone 2011 [56] found matric potential at field capacity ψCC = −6.5 kPa, corresponding to θˆ CC, for τ of 0.01 of KS, in a sandy texture soil in the Cerrado. The moisture at the inflection point (θˆ IP) is considered optimal for soil preparation [29] and correlates well with moisture at field capacity [56]. At the inflection point of the retention curve the first derivative changes from negative (convexity) to positive (concavity), and the second derivative is zero [57], which is relevant to determine the estimation of θˆ CC.

An analysis of the moisture at the inflection point of the soil water retention curve (Table 5; Figure 2) indicated that the drainage rate of the soil corresponded to *p* = 0.050 of the KS. The higher magnitude of the drainage rate of the analyzed soil was due to the texture, predominantly of the coarse sand fraction. For *p* = 0.050, there was ψCC of −4.602 kPa (0 to 0.2 m), −5.711 kPa (0.2 to 0.4 m), and −5.942 kPa (0.4 to 0.6 m). Regarding the time for field capacity to occur (tCC measured in hours), the soil showed rapid vertical drainage due to the magnitude of the saturated hydraulic conductivity (KS). In this case, the high magnitude of K<sup>S</sup> influenced the low values of tCC. The sand fraction of the soil also interferes with the matric potential corresponding to the permanent wilting point.

A higher magnitude of pore diameter was observed in the topsoil layer due to the higher coarse sand content (Table 1). The influence of the coarse sand fraction on the larger pore diameter and, consequently, lower water availability in predominantly sandy soils was also found by Fidalski et al., 2013 [6]. In soils with a predominance of sand, especially coarse sand, soil aggregation is reduced, not allowing the development of a more diverse pore network, especially in the smaller diameter pores, predominantly formed by organomineral interactions in the clay fraction [9,16]. Importantly, larger pore diameters require lower energy for soil water removal [9], consistent with results obtained for saturated hydraulic conductivity Ks and unsaturated hydraulic conductivity K(θ) (Figure 1).

#### 4.1.4. Water Availability to Plants

Fidalski et al., 2013 [6] studied soils of the Paranavaí Sandstone Formation, in the northwestern region of Paraná and found that soil use and management interfered with water availability to plants. Fidalski et al., 2013 [6] also found available water capacity (AWC) in the layer between 0 and 0.4 m of 48.60 mm and 58.70 mm in soil under croplivestock integration and pineapple, respectively. Helbel Junior and Fidalski 2017 [58] found an average water storage capacity of 1 mm/cm in the 0 to 0.6 m layer, in *Latossolo Vermelho distrófico* (Oxisol), medium texture, in the northwestern region of Paraná.

In the 0 to 0.6 m soil layer, the following were verified: AWC<sup>1</sup> = 68.78 mm; AWC<sup>2</sup> = 78.67 mm; EAW<sup>1</sup> = 64.03 mm and EAW<sup>2</sup> = 73.92 mm. AWC<sup>2</sup> overestimates soil water availability because the actual volumetric moisture at the permanent wilting point of the analyzed soil occurs before –1500 kPa (θPMP). EAW<sup>1</sup> underestimates soil water availability because the actual volumetric moisture at field capacity occurs at matric potentials lower than –6 kPa (θCC), highlighting the importance of the actual estimation of volumetric moisture at field capacity [56] and the hydraulic cut-off point [25].

Low water availability to plants in sandy soils can provide the occurrence of water deficit during the sugarcane cycle, as evidenced by Gurski et al., 2020 [59] and Araújo 2019 [60]. High soil hydraulic conductivity (K(θ)) and low time to reach moisture at field capacity (tCC) intensify the occurrence of water deficit. Gurski et al., 2020 [59] verified that water deficiency and surplus were concentrated mainly in development phase II (vegetative growth), in the 1997-98 to 2008-09 crops, in Paranavaí, Paraná, requiring rescue irrigation. Thus, the physical-hydric attributes determined in this study should be considered in the management of irrigation and the plant in each phase of its development.

The effect of sand contents on the water retention of the Caiuá Sandstone, from the northwestern Paraná region, was confirmed with the inverse correlation of coarse sand contents with microporosity, θPMP, clay, fine sand, and θIP, agreeing with [42]. The results again indicated the importance of a more detailed evaluation of the physicalhydric attributes of sandy soils, especially for monitoring the physical quality and water management in the analyzed soil.

#### *4.2. Sustainability Implications*

Soil type and properties at different depths influence the development of sugarcane [61]. While a shallow groundwater table can impact soil moisture dynamics [62], this research was conducted in an area that has a water table that is from 5 to 50 m deep [18], which is associated with vertical, subsurface soil water movement [17]. Yield modeling of sugarcane involves a detailed understanding of the impact of soil water storage dynamics on crop yield in Brazil. This is especially important for soils with low water holding capacity across water deficit gradients which have been modeled in northwestern São Paulo state using a crop water balance model combined with GIS [63]. In Köppen classification Cfa climates like in this study, the average yield response to full irrigation was 22 t/ha using the FAO-AZM simulation model [64]. Greater yield responses were in sandy soils with lower soil water holding capacities, where moisture from rainfall concentrated at the start of the dry season is not retained for as long as other soil types [64]. This research improves understanding of water dynamics for sugarcane grown in such sandy soils and how this relates to crop yield.

More accurate sugarcane models for Brazil's sandy soils are essential to optimize water management decisions to (1) invest in irrigation in order to increase crop productivity and to (2) withhold water, especially in irrigated systems prior to harvest to maximize sugar content. Yield gaps for sugarcane can be closed by using crop breeding for droughttolerant cultivars [65] and irrigation [65,66]. If future climate changes increase sugarcane yield, payback for irrigation investments could decrease from 14 years to only 3 years, making such investments more feasible [67]. However, irrigation investments may be challenging since medium-sized sugarcane farms tend to be the most profitable [68], with potentially less capital available for such upfront purchases. Dias and Sentelhas 2018 [69] simulated sugarcane with APSIM-Sugar, showing stalk biomass reductions of 4% are ideal for maximizing sucrose yields. This can be accomplished by withholding water for drying off prior to harvest for 17 to 31 days for sandy soils with low water holding capacity [69].

The sustainable intensification of sugarcane, which includes increasing yield, has potential environmental benefits in addition to increased crop productivity. Sugarcane has been used recently to rehabilitate degraded pastures in Brazil. Luz et al., 2020 suggest that aside from compaction in the 10 to 20 cm soil layer, conversion from degraded pasture to sugarcane did not result in more soil degradation, whereas conversion from native habitat to extensive pasture resulted in severe soil degradation over time [70]. However, Dias et al., 2021 used APSIM-Sugar to show that climate change from RCP8.5 and RCP4.5 over the next century could result in sugarcane yield reductions in Brazil primarily due to water stress even with irrigation [64]. Future research is needed to improve the accuracy of sugarcane simulation models in sandy soils by validating models developed from experimental results from long-term sugarcane experiments with process-based sugarcane simulation models such as DSSAT-CANEGRO and APSIM-Sugar [71].

#### **5. Conclusions**

In this study, the main physical and hydric attributes of sandy soil were investigated for sugarcane cultivation in southern Brazil. Our results can guide future breeding programs to select sugarcane cultivars more tolerant to water deficit. In sandy soils with a total sand content close to 800 g/kg, the distribution of both coarse sand and fine sand fractions

should be measured to evaluate soil physical and hydric attributes, not just total sand content. The particle size of coarse sand ranges from 2 mm down to 0.2 mm while fine sand ranges from 0.2 to 0.05 mm. Higher coarse sand content increases saturated and unsaturated soil hydraulic conductivity, maximum pore diameter, and macroporosity, as well as reducing microporosity. A higher proportion of fine sand reduces conductivity and increases water retention, especially in deeper soil layers 0.4 to 0.6 m below the soil surface. The actual matric potential at field capacity was −4.60 kPa which is below the matric potential of −6 kPa commonly found in sandy soils in the topmost soil layer. The soil hydraulic cut-off point indicated that the matric potential of the permanent wilting point for sugarcane is reached at −1030 kPa, not −1500 kPa.

Although soil water retention increases with depth, its availability is greatest in the 0 to 0.2 m soil layer, with changes in the actual moisture at field capacity and permanent wilting point. The drainage rate of water at the soil surface is high. Therefore, there is a need for a continuous inflow of water into the soil, whether by precipitation or irrigation in order to ensure adequate soil moisture for the sugarcane crop. Sugarcane yield simulation using a locally validated model demonstrated that adequate crop yields are possible if there is adequate soil water storage for the sandy soils evaluated in this study. Future research can use actual moisture at field capacity and the permanent wilting point to develop a model for estimating sugarcane yield, considering the variation in soil water availability during the growth phase of sugarcane (phase II). This can be used to develop maps of soil water availability and yield estimates in the main sugarcane-producing regions for sandy soils.

**Author Contributions:** Conceptualization, J.L.V. and J.L.M.d.S.; methodology, J.L.V. and A.C.A.; software, J.L.V.; validation, J.L.V.; formal analysis, J.L.V.; investigation, J.L.V.; data curation, J.L.V.; writing—original draft preparation, J.L.V. and J.L.M.d.S.; writing—review and editing, J.L.V.; J.L.M.d.S., A.C.A., R.A.d.O., A.K.H., D.C.d.A., W.M.d.S. and R.M.A.; supervision, J.L.M.d.S.; project administration, R.A.d.O.; resources, A.K.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received funding support from Capes for the Doctoral Program of the lead author.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Study data can be obtained by request to the corresponding author or the first author, via e-mail. It is not available on the website as the research project is still under development.

**Acknowledgments:** We would like to thank RIDESA for the availability of the experimental area and contribution to the research. We would also like to thank the Post-Graduate Program in Soil Science (UFPR), the Modeling of Agricultural Systems Laboratory (LAMOSA)/SCA/UFPR, and C. F. de Araújo-Júnior from the Rural Development Institute (IDR), Londrina, Paraná State. We are also grateful for the comments and edits from three anonymous reviewers which improved the quality of this work.

**Conflicts of Interest:** The authors declare no conflict of interest. Supporting entities had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **Climate Change, Farm Irrigation Facilities, and Agriculture Total Factor Productivity: Evidence from China**

**Hai Li and Hui Liu \***

School of Economics, Hunan Agricultural University, Changsha 410128, China **\*** Correspondence: liuh1220@163.com

**Abstract:** Due to the trend of global warming, individuals from all walks of life have paid close attention to how climate change affects food security. China is a sizable nation with a rich climate and a diverse range of food crops that are of interest to researchers. Additionally, there is little mention of agricultural technology and farm irrigation facilities in academic research on climate change and agricultural economic growth in China. As a result, this study uses the SBM model, panel fixed effect model, and SYS-GMM model to examine the development trend of climate change and food security based on the panel data of Chinese provinces from 2000 to 2020. The study found that China has maintained an average annual growth rate of 4.3% in agricultural total factor productivity (TFP) in recent years, despite the impact of extreme weather. The average annual precipitation has a depressing influence on the TFP in agriculture, while the average annual temperature has the opposite effect. The farm irrigation facilities and agricultural technology's moderating impact is mostly shown in how well they attenuate the impact of climate change on the TFP in agriculture. Food crops have thereby improved their ability to survive natural risks and attain higher yields as a result of advancements in agricultural technology and increasing investment in contemporary farm irrigation facilities. The study's conclusions are used in the article to make the suggestion that strengthening climate change adaptation is necessary to ensure food security. The strategic policy of "storing grain in technology and storing grain in the soil" and the advancement of contemporary agricultural technology must be put into reality while the management system for grain reserves is being improved.

**Keywords:** climate change; farm irrigation facilities; agriculture total factor productivity (TFP); technical advancement

**1. Introduction**

The sustainability of agricultural development, which is essential for human survival, has been severely threatened by climate change [1]. The average annual surface temperature in China is increasing due to global warming at a rate of 0.23 ◦C every ten years. The distribution of temperature and precipitation over space and time will continue to vary due to climate change, which will also increase the frequency and severity of extreme events, including torrential rainstorms, floods, droughts, and insect outbreaks [2]. China has been able to feed 22% of the world's people on only 8% of the planet's territory while using a high-input, high-pollution agricultural growth model, but at a significant cost to resources and the environment. China's agriculture will need to adapt to a more effective, resource-efficient, and environmentally friendly development model in the future against the backdrop of high-quality and sustainable agricultural development as agricultural modernization progresses. Increasing agriculture's total factor productivity (TFP), which is often calculated as the ratio of the total agricultural output to total factor input, is the key to maintaining agricultural economic growth [3]. However, because of how it affects agricultural production and input levels, climate change, particularly extreme weather, has raised a great deal of uncertainty regarding the improvement of the agricultural TFP,

**Citation:** Li, H.; Liu, H. Climate Change, Farm Irrigation Facilities, and Agriculture Total Factor Productivity: Evidence from China. *Sustainability* **2023**, *15*, 2889. https://doi.org/10.3390/ su15042889

Academic Editor: Aaron K. Hoshide

Received: 16 December 2022 Revised: 31 January 2023 Accepted: 2 February 2023 Published: 6 February 2023

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

making it necessary to find a solution [4]. Smallholder farmers can adopt a variety of adaptive strategies in response to climate change, including crop restructuring, more irrigation, and variety selection. However, it is insufficient to rely solely on farmers to take adaptive action [5]. Therefore, it is only possible to make up for the shortcomings of smallholder farmers in adapting to climate change by properly understanding the policy perspective.

Academics have conducted a great deal of research on agricultural TFP, farm irrigation systems, and climate change. First off, there are numerous studies on the effects of climate change on agricultural output in the body of existing literature, and this literature holds a dominant position [6–9]. Additionally, the effects of climate change on agricultural output are mainly detrimental [10]. There is still no agreement among academics regarding the effects of climate change on agricultural productivity due to the significant geographical differences in these effects. In addition, the model that uses unit yield as the explanatory variable allows researchers to examine how climate change affects various crop yields, but it falls short in its understanding of input-output efficiency [11,12].

Based on the difference in climate change indicators, Villavicencio et al. [13] examined the impacts of climate change on the TFP in U.S. agriculture from two perspectives: temperature and precipitation. The results show that annual precipitation had a significant positive effect on the TFP, but the precipitation density had a significant negative effect on the TFP, and temperature change did not have a significant impact on the TFP in most regions. Liang et al. [14] analyzed the effects of climate change on the TFP in agriculture at the seasonal and regional levels, showing that temperature and precipitation in different agricultural regions and seasons accounted for 70% of the TFP changes in U.S. agriculture from 1981 to 2010.

Based on the agricultural TFP of different industries, crops, or cash crops, further research on the effects of climate change on rice yield per plant (TFP) in Japan by Kunimitsu et al. [15] revealed that climate change has different influences on the TFP in different regions. According to Chen and Gong [16], in the short term, extremely high temperatures will negatively affect China's planting industry's TFP, leading to a greater loss in land output. The usual growth cycle of food crops is shortened by a rise in temperature, which also concerns the food supply and consumption due to aberrant precipitation and incalculable harm to grain production per unit area [17,18]. Due to changes in rainfall, evaporation, runoff, and other water-related processes, as well as the subsequent redistribution of water resources over time and geography and subsequent changes in soil moisture, food production is impacted by inadequate water supplies [19,20].

Second, many techniques are now available to measure the change in the agricultural TFP at various phases. The measurement results differ slightly as a result of the different input–output indicators and periods chosen, although the TFP is infrequently considered in studies of farm irrigation facilities. Numerous studies have been conducted on farm irrigation facilities, the majority of which are performance-focused and use the economic growth model as their theoretical framework. These studies examine the contribution of farm irrigation facilities to agricultural output with a focus on the associations between these facilities and economic growth, grain production, farmers' income, and the environment [21]. Scholars build performance evaluation index systems of governance and assess their direct performance, indirect performance, and total performance using the DEA model and network analysis following pertinent evaluation criteria such as the 3E, 4E, and IOO models [22]. The DEA Tobit two-step approach, S-SBM model, Malmquist–Luenberger index, three-stage DEA model, and UHSBM model were derived to provide an empirical analysis of infrastructure supply and investment performance from both static and dynamic perspectives, respectively. These models address the shortcomings of unintended outputs, environmental factors, and random factors.

Additionally, the impacts of farm irrigation infrastructure on agricultural production are mostly seen in four areas: agricultural growth, lowering production costs, fostering structural adjustment in the agricultural business, and fostering agricultural technology advancement [23,24]. The Cobb–Douglas production function model typically includes variables representing infrastructure investment for assessing the marginal effects of farm irrigation facilities on agricultural development [25]. Some academics have recently concentrated on the ecological and production value of agricultural production-related farm irrigation facilities [26]. We can assess the value of ecosystem services in terms of controlling greenhouse gas emissions and controlling the climate as well as conserving water, soil, biodiversity, and the environment to further investigate the effects of farm irrigation facilities on agricultural environmental efficiency [27–29].

The results of the literature search revealed that the majority of academics have also talked about how climate extremes affect food production. Furthermore, it has been established that boosting the building of farm irrigation facilities can significantly raise the quantity and quality of food produced. The relationship between climate change, farm irrigation facilities, and TFP in agriculture, however, has not received much attention from researchers. In addition, many academics have overlooked the impacts of both undesirable outputs and climate change when measuring the TFP in Chinese agriculture. Few academics have concentrated on the effect of farm irrigation facilities on the TFP in agriculture, notably the significance of agricultural technology, when researching the growth processes of China's agricultural economy.

The input structure varies greatly across China's huge territory and wealth of resources, climatic change, and farm irrigation facilities. Additionally, "living off the sky" continues to be the norm in most places. It is clear that two factors need to be prioritized to improve the agricultural TFP. To start, there are differences in the initial factor input for farm irrigation facilities, the structure of endowments, and the rate at which various regions develop agricultural technology. The second is the fluctuation of climatic conditions. The production structure of agricultural economic development changes with the environment. Thus, we need to precisely study how investments in farm irrigation facilities and climate change affect the trajectory of the agricultural TFP.

In light of this, using inter-provincial panel data and a theoretical model of economic growth, this study undertakes an empirical test. The SBM model is utilized in this work to assess the influence of undesirable production on China's agricultural TFP. In addition, the panel fixed-effect model and SYS-GMM model are suggested to precisely investigate the relationship between climate change, farm irrigation facilities, and China's agricultural TFP. On the one hand, we examine the crucial part that farm irrigation facilities play in the process of how harsh climate impacts food production. On the other hand, in order to supplement the already-present data on the variables impacting TFP in agriculture, the crucial component of technological advancement is introduced. It offers practical policy recommendations for managing the shocks brought on by climate change and guaranteeing food security.

This study looks at the effects of climate change, agricultural technology, and farm irrigation facilities on agricultural TFP and aims to innovate in the following areas:

(1) The impact of extreme weather on agricultural production. With global warming, extreme meteorological disaster events such as heavy rainfall, floods as well as droughts are frequent. The annual average surface temperature in China is rising at a rate of 0.26 ◦C/decade, bringing unpredictable changes to the spatial and temporal distribution of temperature and precipitation and increasing challenges to sustainable agricultural development. Therefore, this paper chooses the average annual precipitation and average temperature as the indicators of climate change, which is different from some scholars.

(2) The impact of non-desired output, i.e., agricultural carbon emissions, is taken into account when measuring TFP in agriculture. In addition, this method is different from other studies, such as DEA.

(3) Introducing two intermediary variables—farm irrigation facilities and agricultural technology—to broaden the scope of the research.

(4) Using the most recent years of data from 2000 to 2020 for the agricultural TFP measurement. Furthermore, this paper provides a scientific foundation and policy suggestions for future food security. provinces in China. The study area is shown in Figure 1. There are five different types of terrain in China's general topography: mountains, plateaus, basins, plains, and hills. The country is high in the west and low in the east. This

(4) Using the most recent years of data from 2000 to 2020 for the agricultural TFP measurement. Furthermore, this paper provides a scientific foundation and policy sug-

To this end, this paper is organized as follows: Section 2 illustrates the study area, methodology, models, indicator selection, and data sources used in this study. Section 3 discusses theoretical analysis and analyzes the impact of climate change and farm irrigation facilities on the total factor productivity growth in agriculture using an econometric

The study area of the article is located in the eastern part of Asia and the west coast of the Pacific Ocean. At present, there are 34 provincial administrative regions in China, including 23 provinces, 5 autonomous regions, 4 municipalities directly under the central government, and 2 special administrative regions. Due to the limitation of data availability, this part of the empirical sample data mainly comes from relevant statistics for 31

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 4 of 19

model. Section 4 contains research conclusions and policy implications.

gestions for future food security.

**2. Materials and Methods** 

*2.1. Study Area* 

To this end, this paper is organized as follows: Section 2 illustrates the study area, methodology, models, indicator selection, and data sources used in this study. Section 3 discusses theoretical analysis and analyzes the impact of climate change and farm irrigation facilities on the total factor productivity growth in agriculture using an econometric model. Section 4 contains research conclusions and policy implications. offers a range of possibilities and settings for the growth of China's industry and agriculture. In China, the north and south experience significantly different wintertime temperatures, and summertime highs are typical across the board. Furthermore, the distribution of yearly precipitation is more southerly than northern, with more rain falling in the summer and autumn and less in the winter and spring. Agricultural productivity and other

#### **2. Materials and Methods** things are intimately tied to rainfall conditions. Due to its size and the huge variations in its temperature and precipitation patterns,

#### *2.1. Study Area*

The study area of the article is located in the eastern part of Asia and the west coast of the Pacific Ocean. At present, there are 34 provincial administrative regions in China, including 23 provinces, 5 autonomous regions, 4 municipalities directly under the central government, and 2 special administrative regions. Due to the limitation of data availability, this part of the empirical sample data mainly comes from relevant statistics for 31 provinces in China. The study area is shown in Figure 1. China has a complex and varied climate. This makes China a good location for growing the majority of the world's crops. China's climate has several elements that are helpful for the growth of agricultural production, but it also has some disadvantages. Droughts, floods, cold waves, and typhoons are the most common severe weather phenomena that have a significant influence on China. These disasters frequently have a negative impact on agricultural productivity and farmers' lives.

**Figure 1**. Digital elevation model of the study area. (**a**) Average annual precipitation in China from 2000 to 2020. (**b**) Prevalence of irrigation in China from 2000 to 2020; the ratio of the effective irrigated area to the cultivated area is known as the prevalence of irrigation. The figure includes the provincial boundaries (gray line). The files used to create the map are licensed under GS (2019) 1822 **Figure 1.** Digital elevation model of the study area. (**a**) Average annual precipitation in China from 2000 to 2020. (**b**) Prevalence of irrigation in China from 2000 to 2020; the ratio of the effective irrigated area to the cultivated area is known as the prevalence of irrigation. The figure includes the provincial boundaries (gray line). The files used to create the map are licensed under GS (2019) 1822 available at: http://www.gov.cn (accessed on 1 December 2022). The data on annual precipitation are taken from the National Weather Data Network (http://data.cma.cn, accessed on 1 December 2022). The data on the effective irrigated area to the cultivated area are taken from the EPS database (https://www.epsnet.com.cn, accessed on 1 December 2022), the *China Statistical Yearbook*, and the *China Rural Statistical Yearbook* (http://www.stats.gov.cn, accessed on 1 December 2022).

There are five different types of terrain in China's general topography: mountains, plateaus, basins, plains, and hills. The country is high in the west and low in the east. This offers a range of possibilities and settings for the growth of China's industry and agriculture. In China, the north and south experience significantly different wintertime temperatures, and summertime highs are typical across the board. Furthermore, the distribution of yearly precipitation is more southerly than northern, with more rain falling in the summer and autumn and less in the winter and spring. Agricultural productivity and other things are intimately tied to rainfall conditions.

Due to its size and the huge variations in its temperature and precipitation patterns, China has a complex and varied climate. This makes China a good location for growing the majority of the world's crops. China's climate has several elements that are helpful for the growth of agricultural production, but it also has some disadvantages. Droughts, floods, cold waves, and typhoons are the most common severe weather phenomena that have a significant influence on China. These disasters frequently have a negative impact on agricultural productivity and farmers' lives.

#### *2.2. Methods*

The impact of undesirable outputs is not taken into account when measuring efficiency using the conventional data envelopment analysis (DEA), which only considers economic benefits. This ignores the issue of input–output slackness and is at odds with the actual agricultural production process. This study employs an over-efficient SBM model with non-desired outcomes, which explicitly incorporates the slack variables of each input and output into the objective function. The impact of the slack variables on the measured values is discussed as the total factor productivity (TFP) of each province in the nation is measured using the Max DEA program from 2000 to 2020. The model expressions are as follows:

$$\begin{split} \rho^{\*} &= \min \frac{1 - \frac{1}{m} \sum\_{i=1}^{m} \frac{S\_i^{-}}{z\_{i0}}}{1 + \frac{1}{S\_1 + S\_2} \left(\sum\_{j=1}^{S\_1} \frac{S\_j^{\delta}}{y\_{j0}^{\delta}} + \sum\_{k=1}^{S\_2} \frac{S\_k^{b}}{z\_{k0}^{\delta}}\right)} \\ \text{s.t.} & \begin{cases} \mathbf{x}\_0 = \mathbf{X}\boldsymbol{\lambda} + \mathbf{S}^-, y\_0^{\mathcal{S}} = \mathbf{Y}^{\mathcal{S}}\boldsymbol{\lambda} - \mathbf{S}^{\mathcal{S}}, z\_0^b = \mathbf{Z}^b \boldsymbol{\lambda} + \mathbf{S}^b \\ \mathbf{S}^- \ge \mathbf{0}, \mathbf{S}^{\mathcal{S}} \ge \mathbf{0}, \mathbf{S}^b \ge \mathbf{0}, \boldsymbol{\lambda} \ge \mathbf{0} \end{cases} \end{split} \tag{1}$$

where, in the formula, *ρ* ∗ is the agricultural TFP, where 0 < *ρ* <sup>∗</sup> ≤ 1; *S* −, *S g* , and *S <sup>b</sup>* are the slack vectors of inputs, desired outputs, and undesired outputs, respectively. *x<sup>i</sup>* , *y g j* , and *z b k* are the input of *i*, the desired output of *j*, and the non-desired output vector of *k*, respectively. "0" is the evaluated unit. *m*, *S*1, and *S*<sup>2</sup> are the number of input, desired output, and non-desired output elements. *X*, *Y g* , and *Z <sup>b</sup>* are matrices consisting of inputs, desired outputs, and undesired outputs. *λ* is the weight vector. When *S* − = *S <sup>g</sup>* = *S <sup>b</sup>* = 0, *ρ* ∗ = 1, signifying that the decision unit is totally legitimate; otherwise, it denotes a loss and necessitates adjusting the input and output quantities.

Malmquist's proposal for the Malmquist Index was merged with the DEA theory to create the intertemporally variable TFP in 1953, and the equation is as follows:

$$TFP = \left[\frac{D\_0^{t+1}(\mathbf{x}\_{t+1}, \mathbf{y}\_{t+1})}{D\_0^{t+1}(\mathbf{x}\_t, \mathbf{y}\_t)} \times \frac{D\_0^t(\mathbf{x}\_{t+1}, \mathbf{y}\_{t+1})}{D\_0^t(\mathbf{x}\_t, \mathbf{y}\_t)}\right]^{1/2} = E\mathbf{C} \times T\mathbf{C} \tag{2}$$

where the technical efficiency change index and the technical progress change index, respectively, are denoted by the letters EC and TC. The expressions for each of these two are as follows:

$$\begin{aligned} EC &= SEC \times PEC = \frac{S\_0^t(\mathbf{x}\_t, y\_t)}{S\_0^t(\mathbf{x}\_{t+1}, y\_{t+1})} \times \frac{D\_0^t(\mathbf{x}\_{t+1}, y\_{t+1})}{D\_0^t(\mathbf{x}\_t, y\_t)} \\ TC &= \left[ \frac{D\_0^t(\mathbf{x}\_{t+1}, y\_{t+1})}{D\_0^{t+1}(\mathbf{x}\_{t+1}, y\_{t+1})} \times \frac{D\_0^t(\mathbf{x}\_t, y\_t)}{D\_0^{t+1}(\mathbf{x}\_t, y\_t)} \right]^{1/2} \end{aligned} \tag{3}$$

where, in the formula, *SEC* and *PEC* are the scale efficiency change index and pure technical efficiency change index, respectively. The criteria for determining *EC*, *TC*, *SEC*, and *PEC* are the same. TFP > 1 denotes that the total factor productivity increases, TFP < 1 denotes that the total factor productivity decreases, and TFP = 1 denotes that the total factor productivity does not cause changes.

#### *2.3. Model Specification*

The article analyzes the data by constructing an adjustment model and mediation model using the panel fixed effect method and the SYS-GMM method. The text that follows displays the model's precise form.

#### 2.3.1. Adjustment Model

The agricultural TFP is influenced by climate change and farm irrigation facilities. In order to study the adjustment effect of agricultural technology, the interaction term was set separately for precipitation and temperature using it. And, build the following model.

$$\begin{aligned} TFP\_{i,t} &= \mathfrak{a}\_0 + \mathfrak{a}\_1 pre\_{i,t} + \mathfrak{a}\_2 tem\_{i,t} + \mathfrak{a}\_3 fru\_{i,t} + \mathfrak{a}\_4 \mathfrak{x}\_{i,t} + \mu\_i + \mathfrak{e}\_{i,t} \\ TFP\_{i,t} &= \mathfrak{a}\_0 + \mathfrak{a}\_1 pre\_{i,t} + \mathfrak{a}\_2 tem\_{i,t} + \mathfrak{a}\_3 fru\_{i,t} + \mathfrak{a}\_4 tc\_{i,t} + \mathfrak{a}\_5 tc\_{i,t} \* pre\_{i,t} \\ &+ \mathfrak{a}\_6 tc\_{i,t} \* term\_{i,t} + \mathfrak{a}\_7 \mathfrak{x}\_{i,t} + \mu\_i + \mathfrak{e}\_{i,t} \end{aligned} \tag{4}$$

#### 2.3.2. Mediation Model

A dynamic process affecting the agricultural TFP is influenced by beginning conditions as well as current factors affecting the agricultural TFP, among other things. In order to determine how climate change and farm irrigation facilities affect the agricultural TFP, this study attempts to control the initial conditions, incorporate the lagging term of the agricultural TFP (*TFPi*,*t*−<sup>1</sup> ) into the regression model, and build the following model.

$$\begin{aligned} TFP\_{i,t} &= a\_0 + a\_1 pre\_{i,t} + a\_2tem\_{i,t} + a\_3fru\_{i,t} + a\_4x\_{i,t} + \mu\_i + \varepsilon\_{i,t} \\ TC\_{i,t} &= a\_0 + a\_1pre\_{i,t} + a\_2tem\_{i,t} + a\_3fru\_{i,t} + a\_4x\_{i,t} + \mu\_i + \varepsilon\_{i,t} \\ TFP\_{i,t} &= a\_0 + a\_1pre\_{i,t} + a\_2tem\_{i,t} + a\_3fru\_{i,t} + a\_4tc\_{i,t} + a\_5x\_{i,t} + \mu\_i + \varepsilon\_{i,t} \\ TFP\_{i,t} &= a\_0 + a\_1TFP\_{i,t-1} + a\_2cli\_{i,t} + a\_3term\_{i,t} + a\_4fru\_{i,t} + a\_5tc\_{i,t} + a\_6x\_{i,t} + \mu\_i + \varepsilon\_{i,t} \end{aligned} \tag{5}$$

In Formulas (4) and (5), *TFPit* is the agricultural TFP of province *i* in year *t*; *preit* represents the climate variable of province *i* in year *t*, that is, the yearly precipitation; *temit* is the average temperature; *f ruit* is the input of farm irrigation facilities in year *t* in province *i*, that is, the effective irrigation area per capita; *tcit* is the input variable for the Malmquist–Luenberger index in year *t* in province *i*; *xit* is a control variable, including the rural road density, agricultural structure, population density, and fiscal decentralization; α is the parameter to be estimated; *µ<sup>i</sup>* represents the fixed effect of each province; and *εi*,*<sup>t</sup>* is the random error term. All variables are logged to avoid heteroskedasticity in the model.

#### *2.4. Variables*

#### 2.4.1. The Agricultural TFP

Output variables: Including the expected output and undesired output, the expected output makes use of the grain output, whereas undesired output typically refers to nonpoint-source pollutants such as chemical oxygen demand (COD), nitrogen (N), or phosphorus (P), as well as other pollutants or agricultural carbon emissions. This paper intends to use agricultural carbon emissions as the undesired output.

The sources of agricultural carbon emissions have diverse characteristics. The sources of agricultural carbon emissions have been comprehensively identified as chemical fertilizers, pesticides, agricultural diesel, agricultural film, land plowing, and irrigation power consumption [30,31]. A method for calculating agricultural carbon emissions is built based on the sources of carbon in agriculture that have been recognized:

$$E = \sum E\_{\bar{i}} = \sum T\_{\bar{i}} \times \delta\_{\bar{i}} \tag{6}$$

where, in the formula, *E* is the total amount of agricultural carbon emissions, *i* is the type of agricultural carbon sources, *T<sup>i</sup>* is the consumption of each carbon source, and *δ<sup>i</sup>* is the carbon emission coefficient of each carbon source. The coefficient is 0.8956 kgC/kg, the carbon emission coefficient of pesticides is 4.9341 kgC/kg, the carbon emission coefficient of the agricultural film is 5.18 kgC/kg, the carbon emission coefficient of diesel is 0.5972 kgC/kg, the carbon emission coefficient of tillage is 312.6 kgC/hm<sup>2</sup> , and the carbon emission coefficient of agricultural irrigation should be 25 kgC/hm<sup>2</sup> [32,33].

Input variables: This paper only chooses seven indicators as input variables—labor, machinery, fertilizer, pesticide, agricultural film, diesel oil, and land—in order to comply with the empirical rule of the decision-making unit (DMU) and the number of input variables in DEA analysis [34–36]. The mechanical power input, for example, refers to the total power index of various agricultural machines, including tractors, balers, and seeders used in agricultural output, etc. [37]. The labor force input is based on the number of rural residents per unit area in each region [38]. The amount of agricultural chemical fertilizers applied in each region (in pure volume), pesticide input (the amount of pesticides per unit area in each region), agricultural film input (the amount of agricultural plastic film per unit area in each region), diesel input (the amount of agricultural diesel per unit area in each region), and land input (the effective irrigated area at the end of each region) are used because considering that the article measures the TFP of agriculture, combined with the common phenomenon of agricultural replanting, fallow, and abandonment in China, the ratio of effective irrigated area to cultivated land area and crop sown area are used instead. Agricultural land input is more accurate.

#### 2.4.2. Climate Variables

The influence of climate change on agricultural production is mostly evident in temperature and precipitation [39]. Since this paper focuses on the agricultural TFP of food crop production, annual precipitation is used to measure the impact of climate change on agriculture. The accumulated temperature variable, which is frequently used in agronomy, reflects the impact of temperature on the growth and development of food crops from two aspects: temperature and time [40]. The TFP's impact is comparatively more consistent. In order to quantify the impact of climate change on the agricultural TFP, this research constructs the relationship between annual precipitation and the average temperature and TFP [41].

#### 2.4.3. Farm Irrigation Facilities

At present, scholars usually utilize two types of indicators, monetary and physical, for measurement. However, monetary indicators tend to depart from the true worth of infrastructure since they are usually direct sums of investment volumes, while physical indicators are similarly incorrect due to variances in units of measurement. Usually, the more accurate way is to utilize the perpetual inventory method to estimate agricultural irrigation facilities, but the selection of the depreciation rate and initial capital stock has a significant impact on the inventory results [42,43]. Therefore, experts have varying estimations of farm irrigation facilities. Therefore, based on synthesizing the available research literature, this work selects the effective irrigated area as an index to measure the input of farm irrigation infrastructure.

#### 2.4.4. Control Variables

Several factors affect grain output, farmers' income, and economic development level. In order to objectively evaluate the impact of climate change and farm irrigation facilities on the agricultural TFP, combined with previous research and actual conditions, the model also selected the agricultural structure, and seven disaster factors that have a significant impact on the agricultural TFP are used as control variables (see Table 1).


**Table 1.** Index system for evaluation of agricultural infrastructure governance efficiency.

Note: The table is organized by the author.

#### *2.5. Data Sources and Statistical Analysis*

The time range of the above variables is from 2000 to 2020. The source data are primarily taken from the EPS database (https://www.epsnet.com.cn, accessed on 1 December 2022), the *China Statistical Yearbook*, and the *China Rural Statistical Yearbook* (http: //www.stats.gov.cn, accessed on 1 December 2022).

The economic data of agricultural inputs and outputs as well as climate data for each province from 2000 to 2020 were compiled into a table in accordance with pertinent databases and statistical yearbooks. A descriptive statistical analysis was then carried out for each indicator's data, and the results are shown in Table 2.

**Table 2.** Descriptive statistics of different variables in China from 2000 to 2020.


Note: Data compiled by the authors.

#### **3. Results**

#### *3.1. Trend and Regional Analysis of Agricultural Total Factor Productivity (TFP)*

The agricultural TFP for the entire country is calculated in this research using the Max DEA software and the super-efficiency SBM model, incorporating unexpected production. The dynamic trajectory of the agricultural TFP in various grain-producing regions is considerably diverse from 2000 to 2020, as illustrated in Figure 2. The primary regions for producing grain have an upward tendency overall. The agricultural TFP in this region has been in a good and steady state for a long time, and the agricultural output has resulted in some economic gains, except in 2003, when the governance efficiency was lower than 0.5; the agricultural TFP in the main grain-selling locations varies substantially.

**3. Results** 

*3.1. Trend and Regional Analysis of Agricultural Total Factor Productivity (TFP)* 

The agricultural TFP for the entire country is calculated in this research using the Max DEA software and the super-efficiency SBM model, incorporating unexpected production. The dynamic trajectory of the agricultural TFP in various grain-producing regions is considerably diverse from 2000 to 2020, as illustrated in Figure 2. The primary regions for producing grain have an upward tendency overall. The agricultural TFP in this region has been in a good and steady state for a long time, and the agricultural output has resulted in some economic gains, except in 2003, when the governance efficiency was lower than 0.5; the agricultural TFP in the main grain-selling locations varies substantially.

**Figure 2.** Internal mechanism of farm irrigation facilities affecting the agricultural TFP operation diagram. (China can be divided into three categories, including grain production areas, grain marketing areas, and grain balance areas. Of these, there are 13 grain production areas, including Heilongjiang, Jilin, Liaoning, Inner Mongolia, Hebei, Henan, Shandong, Jiangsu, Anhui, Jiangxi, Hubei, Hunan, and Sichuan; 7 grain marketing areas, including Beijing, Tianjin, Shanghai, Zhejiang, Fujian, Guangdong, and Hainan; and 11 grain balance areas, including Shanxi, Ningxia, Qinghai, Gansu, Tibet, Yunnan, Guizhou, Chongqing, Guangxi, Shaanxi, and Xinjiang). **Figure 2.** Internal mechanism of farm irrigation facilities affecting the agricultural TFP operation diagram. (China can be divided into three categories, including grain production areas, grain marketing areas, and grain balance areas. Of these, there are 13 grain production areas, including Heilongjiang, Jilin, Liaoning, Inner Mongolia, Hebei, Henan, Shandong, Jiangsu, Anhui, Jiangxi, Hubei, Hunan, and Sichuan; 7 grain marketing areas, including Beijing, Tianjin, Shanghai, Zhejiang, Fujian, Guangdong, and Hainan; and 11 grain balance areas, including Shanxi, Ningxia, Qinghai, Gansu, Tibet, Yunnan, Guizhou, Chongqing, Guangxi, Shaanxi, and Xinjiang).

The other 25 Chinese provinces increased positively between 2000 and 2020, as indicated in Table 3, except for Guangdong, Hainan, Guizhou, and Shaanxi, where the TFP was less than 1. The fact that it was larger than 0.970 and that Jilin Province had the highest score and Guizhou Province had the lowest index among them showed that TFP in China was still progressing well. Additionally, the average annual growth rate of the agricultural TFP in the 29 provinces was 4.3%, slightly higher than the findings of the other scholar. The ability to tolerate natural risks has increased grain production in recent years, despite the impact of harsh weather. This is because investments in modern irrigation and water conservation facilities have increased due to technological advancement. However, the precise causes demand further investigation. The other 25 Chinese provinces increased positively between 2000 and 2020, as indicated in Table 3, except for Guangdong, Hainan, Guizhou, and Shaanxi, where the TFP was less than 1. The fact that it was larger than 0.970 and that Jilin Province had the highest score and Guizhou Province had the lowest index among them showed that TFP in China was still progressing well. Additionally, the average annual growth rate of the agricultural TFP in the 29 provinces was 4.3%, slightly higher than the findings of the other scholar. The ability to tolerate natural risks has increased grain production in recent years, despite the impact of harsh weather. This is because investments in modern irrigation and water conservation facilities have increased due to technological advancement. However, the precise causes demand further investigation.

Ningxia and Shanghai had the lowest indexes, both of which were not high, showing that further advancement in agricultural technology is needed to accelerate the development of modern agriculture. Less than half of the provinces reached above 1 on the technological progress and change index (TC) from 2000 to 2020, with the lowest indexes being low values in both cases. Only four provinces had an index of technical efficiency change (EC) below 1, which was consistent with the TFP comparison result and showed that China's agricultural technological level was in good shape. There are 11 provinces in which all three major indexes are larger than 1, and EC is primarily responsible for the rise in TFP, according to the three major indexes. It is clear that in order to raise the local TFP, each province in the nation should concentrate on raising its technological level and increasing its investment in agricultural production. Ningxia and Shanghai had the lowest indexes, both of which were not high, showing that further advancement in agricultural technology is needed to accelerate the development of modern agriculture. Less than half of the provinces reached above 1 on the technological progress and change index (TC) from 2000 to 2020, with the lowest indexes being low values in both cases. Only four provinces had an index of technical efficiency change (EC) below 1, which was consistent with the TFP comparison result and showed that China's agricultural technological level was in good shape. There are 11 provinces in which all three major indexes are larger than 1, and EC is primarily responsible for the rise in TFP, according to the three major indexes. It is clear that in order to raise the local TFP, each province in the nation should concentrate on raising its technological level and increasing its investment in agricultural production.

The national agricultural total factor productivity (TFP) value, which accounts for inputs from farmland and water facilities, can be seen in Figure 3. It fluctuated from 0.821 in 2000 to 1.040 in 2004 and then from 2005 to 2020, except in 2005, 2009, 2010, and 2014, wherein the TFP remained above 1. The average TFP showed an upward trend from 2000 to 2004 and an upward trend from 2005 to 2020, with the smaller fluctuations being brought on by the varying trend of the technological progress change index (TC). The technical efficiency index (EC) is influenced by both the scale efficiency change index (SEC) and the pure technical efficiency change index (PEC), with the mean value of the pure technical efficiency change index (PEC) trending more similarly while the scale efficiency change index has a different trend. Additionally, the mean value of total factor productivity (TFP) and mean value of the technical efficiency index (EC) are relatively similar (SEC). This

suggests that the technical efficiency index (EC) and the pure technical efficiency index of change have a significant impact on the total factor productivity (TFP) in agriculture (PEC).


**Table 3.** TFP of agriculture in China from 2000 to 2020.

Note: Data are calculated by Max DEA software. TFP = EC × TC = PEC × SEC × TC.

#### *3.2. Analysis of the Impact of Climate Change on the Agricultural TFP*

The Hausman test findings for model (1) in Table 4 indicate that the panel fixed effect model is preferable to the mixed regression and random effect models. In terms of the regression results of the explanatory variables, while precipitation has played a negative, but not significant, role on the agricultural TFP, temperature has a significant positive effect on the agricultural TFP. The reason is that precipitation irregularity increases the likelihood of geological disasters such as floods and damage to food production, transportation, and other links, endangering food supply and utilization. The annual precipitation varies greatly across the entire nation and exhibits an upward trend, with four peaks occurring in 2010, 2012, 2016, and 2020, so it is not significant. In comparison to other years, 2016 had much more precipitation, whereas the highest amount of precipitation in 2011 was very low (as shown in Figure 4). Data analysis demonstrates that China's average annual accumulated temperature of 13 to 15 ◦C in the range of floating, in general, is rising. In 2007, 2017, and 2019, there were three relatively obvious annual accumulated temperature peaks, and in 2012, there were the most resources; climate change is the most direct characteristic of climate warming. For crops well north, the winter cold brings beneficial effects, resulting in increased grain production [44].

**Figure 3.** Trend chart of the agricultural total factor productivity (TFP) and its decomposition indicators from 2000 to 2020. (TFP = EC × TC = PEC × SEC × TC, TFP is agricultural total factor productivity, EC is the technical efficiency index, TC is the technological progress change index, PEC is the pure technical efficiency change index, and SEC is the scale efficiency change index). **Figure 3.** Trend chart of the agricultural total factor productivity (TFP) and its decomposition indicators from 2000 to 2020. (TFP = EC × TC = PEC × SEC × TC, TFP is agricultural total factor productivity, EC is the technical efficiency index, TC is the technological progress change index, PEC is the pure technical efficiency change index, and SEC is the scale efficiency change index).

*3.2. Analysis of the Impact of Climate Change on the Agricultural TFP*  The Hausman test findings for model (1) in Table 4 indicate that the panel fixed effect model is preferable to the mixed regression and random effect models. In terms of the regression results of the explanatory variables, while precipitation has played a negative, but not significant, role on the agricultural TFP, temperature has a significant positive effect on the agricultural TFP. The reason is that precipitation irregularity increases the likelihood of geological disasters such as floods and damage to food production, transportation, and other links, endangering food supply and utilization. The annual precipitation varies greatly across the entire nation and exhibits an upward trend, with four peaks occurring in 2010, 2012, 2016, and 2020, so it is not significant. In comparison to other years, 2016 had much more precipitation, whereas the highest amount of precipitation in 2011 was very low (as shown in Figure 4). Data analysis demonstrates that China's average annual accumulated temperature of 13 to 15 °C in the range of floating, in general, is rising. In 2007, 2017, and 2019, there were three relatively obvious annual accumulated temperature peaks, and in 2012, there were the most resources; climate change is the most direct characteristic of climate warming. For crops well north, the winter cold brings beneficial effects, resulting in increased grain production [44]. The moderating factors were decentralized to prevent multicollinearity, and the moderating effect was examined using the panel fixed-effect model. The panel fixed-effect model outperformed the random-effect model, according to the Hausman test. According to model (2)'s findings, the temperature had a positive main effect coefficient on the agricultural TFP, whereas precipitation had a negative main effect coefficient. Agricultural technology served as a moderating factor that attenuated some of the negative effects of rising precipitation on the agricultural TFP. The coefficient of the interaction term between agricultural technology and precipitation was positive but insignificant. With the advancement of science and technology, factors such as breeding technology, greenhouses, and the widespread use of agricultural technology have significantly reduced the impact of natural disasters on food production and improved the security of grain production [45–47]. These factors, along with climate change brought on by a lack of resources such as water, heat, and other resources, have contributed to global warming [48]. The negligible coefficient of the interaction term between agricultural technology and precipitation, however, may have two causes. On the one hand, there is an irregular tendency in the current precipitation situation in China. According to research, a 20% decrease in rainfall will render useless any technology used to boost the productivity of food crops. Agricultural technology has a limited ability to reduce precipitation. On the other hand, the regulation effect of irrigation technology on precipitation instability may be weakened due to the early construction and heavy damage of farm irrigation facilities and the serious abandonment of mechanical wells caused by the beginning of the groundwater protection plan, which is consistent with the small regression coefficient of the interaction term between climate change and the input of farm irrigation facilities mentioned above [49].


**Table 4.** Estimation results of climate change on the agricultural TFP.

Note: \*\*\*, \*\*, and \* denote significance levels of 1%, 5%, and 10%, respectively, with standard errors in parentheses; the same below.

**Figure 4.** Trend chart of the average annual precipitation and the average temperature in China from 2000 to 2020 (http://data.cma.cn, accessed on 1 December 2022). **Figure 4.** Trend chart of the average annual precipitation and the average temperature in China from 2000 to 2020 (http://data.cma.cn, accessed on 1 December 2022).

The moderating factors were decentralized to prevent multicollinearity, and the moderating effect was examined using the panel fixed-effect model. The panel fixed-effect model outperformed the random-effect model, according to the Hausman test. According to model (2)'s findings, the temperature had a positive main effect coefficient on the agricultural TFP, whereas precipitation had a negative main effect coefficient. Agricultural

rising precipitation on the agricultural TFP. The coefficient of the interaction term between agricultural technology and precipitation was positive but insignificant. With the advancement of science and technology, factors such as breeding technology, greenhouses, and the widespread use of agricultural technology have significantly reduced the impact of natural disasters on food production and improved the security of grain production [45–47]. These factors, along with climate change brought on by a lack of resources such as water, heat, and other resources, have contributed to global warming [48]. The negligible coefficient of the interaction term between agricultural technology and precipitation, however, may have two causes. On the one hand, there is an irregular tendency in the current precipitation situation in China. According to research, a 20% decrease in rainfall will render useless any technology used to boost the productivity of food crops. Agricultural technology has a limited ability to reduce precipitation. On the other hand, the regulation effect of irrigation technology on precipitation instability may be weakened due to the early construction and heavy damage of farm irrigation facilities and the serious abandonment of mechanical wells caused by the beginning of the groundwater protection plan, which is consistent with the small regression coefficient of the interaction term between climate change and the input of farm irrigation facilities mentioned above [49].

**Table 4.** Estimation results of climate change on the agricultural TFP.

**Variables (1) (2)** 

*lnpre* −0.0160 \*¥(−1.80) −0.00997¥(−1.17) *tem* 0.0264 \*\*¥(2.16) 0.0289 \*\*¥(2.52) *lnfru* 0.108 \*\*\*¥(3.05) 0.0818 \*\*¥(2.47) *lncap* 0.0331¥(1.33) 0.0452 \*¥(1.93) *lnaff* −0.370 \*\*\* −0.266 \*\*\*

*lneng* 0.0563¥(1.10) 0.0275¥(0.57) *lnfer* −0.0756 \*\*¥(−2.06) −0.0661 \*¥(−1.92) *lninv* −0.0177 \*¥(−1.68) −0.0124¥(−1.25) *lnfil* −0.0115¥(−1.10) −0.0143¥(−1.44)

(−6.47) (−4.74)

Agricultural technology, as a moderator variable, significantly reduces the temperature rise's positive influence on the agricultural TFP due to the sharp rise in temperature and aggravation of drought, plant diseases, and insect pests; it also affects the crop quality, causes an increase in extreme weather, and poses a threat to agricultural production [50]. Agricultural technology and a temperature interaction coefficient under the 1% level are significantly negative.

The paper offers a simple theoretical model to explain the influence of climate change on the agricultural TFP and the function of farm irrigation facilities based on the aforementioned literature review. Agriculture is considered to conform to the standard C-D production function form:

$$\mathbf{Y} = \mathbf{A} \mathbf{F} \text{ (K, L)}\tag{7}$$

where Y stands for the agricultural output; A stands for the agricultural TFP; and K and L stand for the capital and labor inputs for agricultural production, respectively.

In addition, A = *θ*T, that is, A is a function of factor allocation (*θ*) and agricultural technological progress (T). Farm irrigation facilities do, however, have a ceiling effect (*θ* ≤ 1) that reduces the efficiency with which agricultural output elements are distributed. Therefore, T' (f) = 0, *θ* 0 (f) ≥ 0, *θ* <sup>00</sup>(f) ≤ 0, and f stands for farm irrigation facilities.

This paper analyses the impact of climate factors on the agricultural TFP, which is denoted by *cli*. The formula is as follows:

$$\frac{dA}{dcli} = \frac{\partial A}{\partial cli} + \frac{\partial A}{\partial \theta(f)T(f)} \cdot \frac{\partial \theta(f)}{\partial cli} \cdot \frac{\partial T(f)}{\partial cli} = \frac{\partial A}{\partial cli} + \frac{\partial A}{\partial f} \cdot \frac{\partial f}{\partial cli'} \tag{8}$$

The influence of climate change on the agricultural TFP is dependent on two factors: the direct effect of climate on the TFP and the indirect effect of climate change on the TFP, which is the input of farm irrigation facilities on the effectiveness of factor allocation θ and technical advancement.

Conclusion 1: Technological progress plays a regulatory role.

#### *3.3. Analysis of the Impact of Farm Irrigation Facilities on the Agricultural TFP*

The panel fixed effect model is used in this study to investigate the relationship between farm irrigation facilities, technological advancement, and the agricultural TFP based on the intermediary effect. The first-order lagged TFP is used as the instrumental variable in the systematic GMM estimate approach to show how climate change affects TFP, and the data passed the robustness test. Table 5 displays the findings, which demonstrate that model (6) passed the test for instrumental variables. The residual terms only possessed a first-order serial correlation, according to the findings of the AR (1) and AR (2) tests, and there was no second-order autocorrelation. All of the Hansen statistics' P-values were higher than 0.1, proving the validity of the instrumental variables.

According to the results of model (5)'s estimation, the agricultural TFP with one lag period passed the positive significance test at the 10% level in all models, suggesting that capital accumulation in the early stages may not have a positive impact on agricultural economic growth in the later stages and that there may be a phenomenon known as diminishing marginal utility. It has been demonstrated, however, that TFP does exhibit "inertia" in time series. A continuous accumulation adjustment process is being used to improve the agricultural TFP.

The coefficient of farm irrigation facilities in model (3) is significantly positive based on the regression results, indicating that these facilities have a positive spillover effect on grain production growth and that accelerating infrastructure investment is a key strategy for enhancing the agricultural TFP. The accumulation of farm irrigation facilities is the primary internal component of technical advancement, which is consistent with model (4)'s finding that the coefficient of farm irrigation facilities is significantly positive. Model (5) shows that agricultural technology and farm irrigation facilities have a significant positive impact on the agricultural TFP, suggesting that technological advancement has a mediating

effect. Model (6) shows that accelerated technical innovation and progress can boost the improvement of the agricultural TFP. These results are consistent with the estimated results of the fixed effects.


**Table 5.** Estimation results of farm irrigation facilities on the agricultural TFP.

Note: L.TFP is the agricultural TFP of the lag period. \*\*\*, \*\*, and \* denote significance levels of 1%, 5%, and 10%, respectively, with standard errors in parentheses.

This research also investigates the underlying mechanisms through which the investment in farm irrigation facilities impacts the agricultural TFP. To find the solution to this question, we can simultaneously derive the input of agricultural irrigation facilities on both sides of Equation (8):

$$\frac{d^2A}{dclidf} = \frac{\partial^2A}{\partial cli \partial f} + \frac{\partial^2A}{\partial^2 f} \cdot \frac{\partial f}{\partial cli'}\tag{9}$$

The left side of Equation (9), which can be broken down into direct impacts and indirect effects of farm irrigation facility investment, shows how farm irrigation facility investments are used to deal with the consequences of climate change on the agricultural TFP.

The technical effects of farm irrigation facilities primarily consist of two aspects, as depicted in the figure: (1) horizontal effect: assuming that the agricultural technical conditions of T1 do not change, the increase in the stock level of farmland water conservancy facilities (as shown in Figure 5, from f1 to f2) will increase the allocation efficiency of

agricultural production factors, and the agricultural TFP will also gradually rise; (2) growth effect: based on the assumption that exogenous agricultural technical circumstances exist, the agricultural TFP will steadily rise from T1 to T2, while the marginal utility will fall. This means that a higher level of farmland and farm irrigation facilities will correlate to a lower ricultural production factors, and the agricultural TFP will also gradually rise; (2) growth effect: based on the assumption that exogenous agricultural technical circumstances exist, the agricultural TFP will steadily rise from T1 to T2, while the marginal utility will fall. This means that a higher level of farmland and farm irrigation facilities will correlate to a lower ′ and A′. and A0 .*Sustainability* **2023**, *15*, x FOR PEER REVIEW 15 of 19 This means that a higher level of farmland and farm irrigation facilities will correlate to a lower ′ and A′. Conclusion 2: Technological progress plays an intermediary role in the process of TFP affected by farm irrigation facilities.

Conclusion 2: Technological progress plays an intermediary role in the process of

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 15 of 20

5%, and 10%, respectively, with standard errors in parentheses.

prove the agricultural TFP.

mated results of the fixed effects.

both sides of Equation (8):

both sides of Equation (8):

TFP.

TFP.

*AR(1)*−2.84 *AR(1) p-value*0.005 *AR(2)*−0.14 *AR(2) p-value*0.888 *Sargan-test*23.03 *Sargan-test p-value*0.113 Note: L.TFP is the agricultural TFP of the lag period. \*\*\*, \*\*, and \* denote significance levels of 1%,

According to the results of model (5)'s estimation, the agricultural TFP with one lag period passed the positive significance test at the 10% level in all models, suggesting that capital accumulation in the early stages may not have a positive impact on agricultural economic growth in the later stages and that there may be a phenomenon known as diminishing marginal utility. It has been demonstrated, however, that TFP does exhibit "inertia" in time series. A continuous accumulation adjustment process is being used to im-

The coefficient of farm irrigation facilities in model (3) is significantly positive based on the regression results, indicating that these facilities have a positive spillover effect on grain production growth and that accelerating infrastructure investment is a key strategy for enhancing the agricultural TFP. The accumulation of farm irrigation facilities is the primary internal component of technical advancement, which is consistent with model (4)'s finding that the coefficient of farm irrigation facilities is significantly positive. Model (5) shows that agricultural technology and farm irrigation facilities have a significant positive impact on the agricultural TFP, suggesting that technological advancement has a mediating effect. Model (6) shows that accelerated technical innovation and progress can boost the improvement of the agricultural TFP. These results are consistent with the esti-

question, we can simultaneously derive the input of agricultural irrigation facilities on

డడ <sup>+</sup> <sup>డ</sup>మ

The left side of Equation (9), which can be broken down into direct impacts and in-

The technical effects of farm irrigation facilities primarily consist of two aspects, as

direct effects of farm irrigation facility investment, shows how farm irrigation facility investments are used to deal with the consequences of climate change on the agricultural

depicted in the figure: (1) horizontal effect: assuming that the agricultural technical conditions of T1 do not change, the increase in the stock level of farmland water conservancy facilities (as shown in Figure 5, from f1 to f2) will increase the allocation efficiency of agricultural production factors, and the agricultural TFP will also gradually rise; (2) growth effect: based on the assumption that exogenous agricultural technical circumstances exist, the agricultural TFP will steadily rise from T1 to T2, while the marginal utility will fall.

This research also investigates the underlying mechanisms through which the investment in farm irrigation facilities impacts the agricultural TFP. To find the solution to this question, we can simultaneously derive the input of agricultural irrigation facilities on

ௗమ

ௗௗ <sup>=</sup> <sup>డ</sup>మ

డడ <sup>+</sup> <sup>డ</sup>మ

The left side of Equation (9), which can be broken down into direct impacts and indirect effects of farm irrigation facility investment, shows how farm irrigation facility investments are used to deal with the consequences of climate change on the agricultural

The technical effects of farm irrigation facilities primarily consist of two aspects, as depicted in the figure: (1) horizontal effect: assuming that the agricultural technical conditions of T1 do not change, the increase in the stock level of farmland water conservancy facilities (as shown in Figure 5, from f1 to f2) will increase the allocation efficiency of ag-

<sup>డ</sup>మ <sup>∙</sup> డ

డ, (9)

<sup>డ</sup>మ <sup>∙</sup> డ

డ, (9)

ௗమ ௗௗ <sup>=</sup> <sup>డ</sup>మ

**Figure 5.** Diagram of internal mechanisms of farm irrigation facilities affecting agricultural TFP op-**Figure 5.** Diagram of internal mechanisms of farm irrigation facilities affecting agricultural TFP operations.

erations. *3.4. Limitations and Implications*  Conclusion 2: Technological progress plays an intermediary role in the process of TFP affected by farm irrigation facilities.

#### *3.4. Limitations and Implications*

Climate change is one of the shifts the world is going through. According to this paper's analysis of climate change in China from 2000 to 2020, which is also in line with the conclusions of Li et al. [51], China is currently undergoing extreme warming. Alexander et al. [52] conducted research on extreme climate change, and the results showed that 70% of the world's land area exhibits a growing trend toward a continuous decline in the number of cold night days and a continuous increase in the number of warm night days. As a result, the Chinese region fits with the trend of global climate change. The results of earlier Climate change is one of the shifts the world is going through. According to this paper's analysis of climate change in China from 2000 to 2020, which is also in line with the conclusions of Li et al. [51], China is currently undergoing extreme warming. Alexander et al. [52] conducted research on extreme climate change, and the results showed that 70% of the world's land area exhibits a growing trend toward a continuous decline in the number of cold night days and a continuous increase in the number of warm night days. As a result, the Chinese region fits with the trend of global climate change. The results of earlier investigations are more compatible with the finding that the intensity of precipitation increases with some variations in precipitation variability [53].

investigations are more compatible with the finding that the intensity of precipitation increases with some variations in precipitation variability [53]. Precipitation was shown to be the principal growth factor impacting grain output and to have a suppressive effect on grain production when the effects of temperature and precipitation on the TFP in agriculture were compared, which was in line with the findings of Yang et al. [54]. The data analysis results revealed that, in contrast to the findings of other studies, the rise in temperature was accompanied by an increase in the agricultural TFP. Combining the research findings of Liu et al. and Haider et al. [55,56] on the Precipitation was shown to be the principal growth factor impacting grain output and to have a suppressive effect on grain production when the effects of temperature and precipitation on the TFP in agriculture were compared, which was in line with the findings of Yang et al. [54]. The data analysis results revealed that, in contrast to the findings of other studies, the rise in temperature was accompanied by an increase in the agricultural TFP. Combining the research findings of Liu et al. and Haider et al. [55,56] on the impact of temperature change on grain production in a particular country, it is possible to conclude that, despite China's ongoing warming trend, the country's warming climate has advanced the start of the warmer season in spring and postponed the start of the cooler season in autumn. Crops have more time to grow, absorb sunlight, photosynthesize, and change matter as a result. This theoretically increases their production capacity and may, in part, result in larger grain yields.

impact of temperature change on grain production in a particular country, it is possible to To sum up, the analysis of this work still has numerous limitations.

Firstly, there are other elements besides climate change that have an impact on the TFP in agriculture. Human factors are also crucial. In theory, investing in farm irrigation facilities can increase the TFP [57]. The policymakers of farm irrigation facility investment are unable to determine the degree of farm irrigation facility investment that is most

appropriate for all types of climate change challenges in all places due to factors such as food security [58,59]. How to open the "last mile" of agricultural technology promotion, which is related to the development of digital platforms for agricultural technology services, as well as how to allow the most cutting-edge and advanced agricultural science and technology achievements to permeate the village, home, and field, is critical. It is essential to obtain local government spending at this time. As a result, future research will continue to focus more on the processes of harsh climate and the effects of anthropogenic variables on the food supply.

Second, this study measures the effect of climatic extremes on agricultural production using average temperature and precipitation density rather than measuring climate extremes directly. According to preliminary findings, temperature change has a large positive impact on the agricultural output increase, while extreme precipitation weather has a negative impact. The measurement and evaluation of extreme weather events as well as the processes by which they affect the productivity of all factors in agriculture, however, require more in-depth study.

Third, due to the size of the research region, there are also significant variations in hydrothermal conditions, particularly for the many food crops that have distinct spatial distribution patterns. Regions or types of changes affected directly by climate change vary greatly [60]. For instance, whereas precipitation increases in South China may increase the likelihood of agricultural disasters, precipitation increases in Northwest China will greatly enhance agricultural output [61,62]. Therefore, future research will concentrate on how extreme climatic change affects the agricultural TFP in various regions of China and for various food crops.

#### **4. Conclusions**

In this study, we examined the mechanism underlying the impact of climate change on the agricultural total factor productivity (TFP), developed a model of the moderating effect, examined the moderating impact of farm irrigation facilities and agricultural technology, and empirically tested the model using panel data from all provinces. The result is the same as that used by Li et al., Alexander et al., and Yang et al. The present paper confirms their findings using the most recent data, different methods, and more scientific indicators of agricultural TFP measurement. The main conclusions are as follows:

(1) Similar to the findings of other studies, this paper finds that climate extremes and farm water facilities have an impact on food output. Agriculture's technological advancements also have a moderating and mediating effect;

(2) In contrast to the findings of other studies, the paper discovers that China's average agricultural TFP between 2000 and 2020 is 4.3%, with a TFP above 1 in 25 provinces, at a somewhat faster rate than other research. Furthermore, this paper makes the case that the average annual cumulative temperature in China varies between 13 and 15 degrees with a general upward trend, and the warming is beneficial for crop growth and preventing freezing calamities, which has a large and positive impact on the agricultural TFP.

The article does, however, have certain shortcomings. Instead of relying solely on data that are available to the public, we must do additional field research for micro-subjects in future studies. Additionally, we must broaden our choice of indicators, particularly for extreme climate change.

The aforementioned findings have significant policy ramifications for combating climate change and guaranteeing domestic food security: (1) in order to handle the food crisis in the short term while strengthening the ability to foresee meteorological disasters, the method for managing grain reserves should be strengthened in light of the constricting effect of climate change on the agricultural TFP; (2) we must aim at the moderating influence of agricultural technology and farm irrigation systems in the process of climate change on the agricultural TFP. On the one hand, we must firmly implement the "storing grain in technology, storing grain in the land" strategy, "short board" swallow the infrastructure, adjust measures to local conditions to adjust the structure of grain cropping systems and

planting, and achieve cultivation optimization of new varieties of crops, taking the comprehensive technology to improve the ability to withstand natural disasters and improve crops' adaptability to environmental changes. On the other hand, although agricultural technology has controlled the precipitation process, its influence has not been substantial, indicating that China's response to climate change and the overall rate of the agricultural technology extension mechanism transformation still need to be improved.

**Author Contributions:** H.L. (Hai Li) generated the conceptualization, methodology, software, and validation; carried out formal analysis, investigation, resources, and data curation; and wrote the original draft, including preparation, review, visualization, supervision, and editing. H.L. (Hui Liu) provided statistical assistance and read, revised, and shaped the manuscript to the present form. H.L. (Hai Li) and H.L. (Hui Liu) co-provided project administration and funding acquisition. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Key Science Fund Project of Hunan Provincial Department of Education, grant number 18A085; Hunan Provincial Philosophy and Social Science Fund Project, grant number 21YBA079; and Hunan Postgraduate Scientific Research Innovation Project, grant number CX2018B403.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** This is not applicable to this article since no datasets were generated.

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

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
