Predicting the influence of extreme temperatures on grain production in the Middle-Lower Yangtze Plains using a spatially-aware deep learning model

Defining Grain Crops.

According to National Bureau of Statistics of P. R. China (2023), the administrative definition of grain crops include rice, wheat, corn, millet, sorghum, and “other grains” (such as barley, oats, buckwheat, etc.). We used grain data from official statistical yearbooks conforming to this definition.

Area of study

Based on examples of previous studies (such as Ding et al. (2017)), the border of the geographical Middle-Lower reaches of the Yangtze River Basin (YRB), and secondary (province-level) administrative boundaries, the Middle-Lower Yangtze Provinces focused in this study are Jiangsu (Administrative code CN032), Zhejiang (CN033), Anhui (CN034), Jiangxi (CN036), Hubei (CN042), and Hunan (CN043) (in Fig. 1). This defined MLYP provincial boundary would also constitute the region from which climate data will be extracted. Shanghai (CN031) is left out of the study not only because its jurisdictional area is a mega-city (larger than municipal-level cities but significantly smaller than provinces) but also due to its apparent inclination toward a financial-centric, non-agricultural economy.

Observation data

Referencing to the typical grain production function proposed by Just & Pope (1979), this study included three provincial-level variables: annual sown area of grain crops(x₁: SownArea) number of primary-industry-employed personnel (x₂: PrimInd, a measure of labor), and total power of agricultural machinery(x₃: MechPower, a measure of physical capital). The number of years elapsed since the first recorded data point implemented in the study (x₄: dt, a time dummy) is also included. dt improves the completeness of agricultural production functions economically as an input elasticities and factor productivity multiplier (Gong, 2018). It also boosts model performance in regression analyses (Just & Pope, 1979).

As shown in Table 1, SownArea, PrimInd, and MechPower are obtained from provincial and China Statistical Yearbooks produced by the National Bureau of Statistics of P. R. China (2022) and cross-checked with Eastern China Statistical Yearbooks from the Eastern China Statistical Information Network (1997). Recent editions of the said Yearbooks are obtained from the respective authorities. The Nanjing Library provides archived historical data.

Four variables are used as climate/meteorological inputs: mean air temperature (tas), measured in °C; total precipitation (pr), measured in mm; HW frequency (hw_freq); and CW frequency (cw_freq). As listed in Table 1B, the GHCN-CAMS Temperature dataset provided by NOAA PSL is used for tas, and pr data is obtained from the CPC Global Unified Gauge-Based Analysis of Daily Precipitation. “CPC Global Unified Temperature” supplies the daily minimum and maximum temperature values used to deduce hw_freq and cw_freq. To maintain cohesion on the temporal axis, the pr dataset is converted to monthly by summing each natural month. On the space axis, climate variables remain as grids.

GCMs

The SSPs outline several different pathways based on varying degrees of socio-economic challenges to mitigation and adaptation. Specifics for each SSP scenario are listed in Table 2. These scenarios range from sustainable development to regional rivalry and fragmentation (Bond & Scott, 2022). SSPs are integral to the ScenarioMIP activity of CMIP6, providing narratives that guide the projection of future emissions and land use in the models.

CMIP6 models and projections used in this study are NASA’s NEX-GDDP CMIP6 dataset prepared by Thrasher et al. (2022). As part of the NASA Earth Exchange (NEX) project, it provides downscaled climate projections from the CMIP6 GCMs across all four “Tier 1” SSPs, making them suitable for regional climate studies. CMIP6 historical and ScenarioMIP experiments SSP126, SSP245, SSP370, and SSP585 of 25 Downscaled GCM models are selected for projecting future climate-induced grain output variability in the MLYP. This dataset is obtained from NASA’s NCCS THREDDS Data Server (NASA NCCS, 2024). More information about these GCMs can be found in Table S5.

Data preprocessing

The proposed model combines spatiotemporal information in reanalyzed meteorological datasets with time series of provincial economic variables to model MLYP grain output. Cases used to train the final model are grouped by year. For each case, gridded meteorological data is first normalized to a uniform resolution, then annual HW and CW frequencies are obtained. To deduce HW and CW frequencies, previously accepted definitions of ETEs, such as that featured in Perkins-Kirkpatrick et al. (2017) and Smith, Zaitchik & Gohlke (2012), are assessed, and, in this study, the 90th percentile of the cumulative density function (CDF) for each point on the 0.5° * 0.5° daily maximum temperature (tasmax) grid 1981-2010 would define the minimum HW-qualifying temperature of a specific location and a day will be labeled as ‘warm day’ when daily tasmax exceeds the threshold over at least three consecutive days. CW is determined following similar principles: a series of cold days on a point in the grid would be considered a coldwave if such a series contains strictly three or more consecutive days with tasmin lower than the 10th percentile value of the CDF for the point’s daily tasmin 1981–2010 data. Specifics regarding ETE definitions can be found in the author’s previous work (Mu & Huang, 2023).

The processed meteorological data are then encoded into eight codes with two distinct ConvAE architectures, each specialized for a specific temporal resolution. Meteorological data are then combined with the corresponding dataset for each economic variable normalized to the closed range [0,1]. Hence, right before the RF regression is performed, each case in the input dataset would contain four economic indices and 32 meteorological indices.

Convolutional autoencoders

Previous studies extensively used convolution, either in the form of a convolutional autoencoder (ConvAE) or a convolutional neural network, to facilitate dimension reduction. Hameed et al. (2022) applied a ConvAE-based deep learning approach for aerosol emission detection. Yu (2022) introduced a constrained dense ConvAE and deep-neural-network-based semi-supervised method to facilitate geological formation recognition. In this study, ConvAEs simplify the final RF regression and improve model performance by reducing potential noise.

Two ConvAE architectures are proposed (Fig. 2, Fig. S2). The spatial convolutional autoencoder (hereafter referred as the sAE) specializes in encoding 2d annual hw_freq and cw_freq datasets. Through two 2d convolution layers (each followed by a Rectified Linear Unit (ReLU) for zero activation) and dense layers, eight-dimension vectors are created as embeddings for each dataset. The spatio-temporal convolutional autoencoder (hereafter referred as the stAE) specialized in encoding 3D monthly tas and pr datasets. Through 3 3D convolution layers (also followed by a ReLU) and dense layers, eight-dimension vectors are created as embeddings for each dataset.

ConvAEs are trained in batches with an unsupervised algorithm based on how well the decoder can reconstruct the sample from the code encoded by the encoder. Mean squared error (MSE) is used to quantify the reconstruction error (i.e., act as the model performance criterion); a higher total MSE loss of the batch penalizes the current weights. Backpropagation computes the gradients. The Adaptive Moment Estimation (Adam) optimizer then uses the computed gradients to update model parameters. For more information on ConvAE training, please reference Supplementary Section S1.1.

Random forest regression

This study utilized the RF regression implementation created by Pedregosa et al. (2011). In RF regression tasks, each tree is constructed by recursively splitting the data based on feature values. Splits are chosen to minimize the variance within each node. Once all trees are grown, the RF regressor aggregates their predictions to form the final model output. In this study, this aggregation is the average of the predictions from all trees. This averaging mitigates overfitting and reduces model variance, leveraging the central limit theorem.

In this study, the RF input set (processed observed economic and climate data in the timeframe [1990,2021]) is split 80%-to-20% into training/validation and testing sets. A grid search optimization process is utilized to tune the four tree-growing hyperparameters of the RF regressor (n_estimators, max_depth, min_samples_split, min_samples_leaf). Additionally, one tree-pruning parameter, ccp_alpha, which quantifies the complexity parameter of the minimal cost-complexity pruning algorithm, is also optimized.

Analyzing the impact of future climate change

Though SSPs contain data regarding land use and agricultural/industrial development, the resolutions of these datasets are not nearly high enough to enable reliable regional studies: GDP and employment projections are included at a national level (Bond & Scott, 2022), and the critical capital factor, power of agricultural machines, is not included in the scenarios. Hence, this study will attempt to evaluate future production possibilities under the current (2021) economic situation, and isolate the effect of climate change for analysis.

Model benchmarking

Projected ETE in MLYP

Before modeling, all climate datasets are pre-processed spatio-temporally. To utilize ScenarioMIP projections of future tasmin and tasmax to deduce ETE frequencies, the thresholds are downscaled linearly to a 0.25° * 0.25° grid, matching the spatial resolution of the NEX-GDDP-CMIP6 GCMs. Then, HW frequencies and CW frequencies are deduced according to the definitions outlined in ‘Data Preprocessing’.

Spatial patterns

Figure 3 shows the spatial distribution of ETEs across three 20-year periods and 2 SSP scenarios. The two extremes differ spatially: HW frequency increase is more significant toward the equator and insignificant in higher latitudes, while CW frequency consistently declines across space–time. Currently, as in the upper left panel of Fig. 3, higher mean annual HW frequencies are observed in the North of the Yangtze river than in the South. Mean HW frequency predicted by selected GCMs for the provinces to the Northern side of the MLYP shows (1) a less significant increase of HW frequency over time and (2) less experiment difference between different SSP scenarios. Among the MLYP provinces, over the last 20 years, Anhui and Zhejiang generally have the highest number of mean HW_freq. However, Anhui’s mean annual heatwave frequency may slightly decrease in the late 21st century, as reflected by heatwave frequency’s Pearson correlation with time graphed in the upper right panels of Fig. 3.

Some discontinuities from observed HW frequency are observed across the spatial distribution of HW frequencies in SSP scenarios. HW hotspots in Northern and Eastern Zhejiang provinces and spots of reduced HW frequencies surrounding major reservoirs in Hunan, Anhui, and Hubei are seen across a few of the 25 GCMs. There is also slight projected decrease in the frequencies of HW occurrence in historically warm regions around MLYP, especially in coastal Zhejiang and Fujian (out of scope).

Temporal patterns

Processed terrestrial HW and CW frequencies are averaged to 4-year means across the MLYP provinces. SSP scenarios all project a general increase in spatial mean HW frequency in Fig. 4A. But, starting with the 2057-2060 group, the SSP370 ensemble is beginning to project a slightly higher number of median HW frequencies than SSP585 (while both still fall under a similar range). This surpassing will be more significant at the end of the century: both the median and maximum of SSP370 ensemble heatwave projections would be approximately 10% more than that of SSP585. This is because of the projected land use changes innate in SSPs (Bond & Scott, 2022). Another observation: the SSP126 scenario projects a much higher level of CW frequency in 2061-2100 than SSP245 (max median difference reaching 2 at 2097-2100, max max difference reaching 1 at 2061-2064). Although 4-year medians and ranges of HW frequency of their groups are relatively similar throughout the projected timeframe, with a numerical max maximum difference merely reaching +1, which did not come until 2097-2100, max of SSP126 annual mean HW frequency projection has also exceeded that of SSP245 in several year groups.

According to Fig. 4B, CW frequency has a stronger negative relationship with time than HW frequencies, especially in scenarios with higher radiative forcings. In SSP585, by the end of the century, more than 50% of the points projects less than two CW events per year. Figure 3 supports this trend. The extent of such decrease does not vary significantly spatially, but there is a clear trend across time. The higher increase in temperature projected by SSP585 made its mean projected CW frequency significantly lower, and Pearson correlation results of SSP585′s projected CW frequency over time have an especially large magnitude that almost doubles that of SSP245′s. Noteworthy temporal trends appear in HW frequency’s temporal projections.

Current grain cultivation in MLYP

In the past decades, Chinese agriculture experienced a transformation from labor-intensive agriculture to industrial agriculture that relies more on machines. As in Fig. 5, this shift from labor- to machine-based agriculture is well visible in the base data implemented in this study: while a high number of personnel are employed in the primary industry in the last decades of the 20th century, none of the MLYP provinces have more than two million workers in agriculture after 2005. Similarly, an increase in the wealth of rural households and a boom of industrialized agriculture contributed to a significant rise in agricultural machinery (Xie & Jin, 2015): total power of agricultural machinery more than tripled in most MLYP provinces over 32 years (1990–2021).

Meteorological data are passed through ConvAEs, and the outputs are compiled together with the four economic variables into 36-dim vectors to form the input dataset of the RF regression. The RF’s parameters are optimized in a search grid (Fig. 6, refer to Section S1.1 for more info). Table 3 shows the importance of each variable in the response. SownArea plays the most crucial role in determining total grain output. This is followed by tas. cw_freq (CW frequency) comes third. The aggregated influence of the eight hwave codes and the eight cwave codes are significantly higher than that of PrimInd, which has a similar level of importance as dt and pr. Observed and predicted total grain production for each province over time are included in Fig. S4.

Model evaluation

The three model evaluation statistics are calculated for both the proposed ConvAE-RF and benchmark models B1-FGLS and B2-RF (Table 4). All three evaluation metrics are used for comparison with the benchmark models. The proposed model has the lowest MAPE (49.2% of B1-FGLS and 59.9% of B2-RF), the highest amount of variance captured (112.0% of B1-FGLS and 105.0% of B2-RF), and the lowest D_γ (19.5% of B1-FGLS and 35.0% of B2-RF, therefore a higher model likelihood) among the three models. Based on the three metrics, the optimized ConvAE-RF model demonstrates superior performance over the multilinear grain production model and the common spatial Random Forest regression model in the MLYP. The second best among the three evaluated models is the common temporal RF regression model. Both RF-based models performed significantly better than the multilinear baseline (B1-FGLS), which aligns with previous research. The results are listed side-by-side in Table S3.

Future grain production in MLYP

Economic data are kept constant at the current level. NEX-GDDP-CMIP6 GCM data are encoded via optimized ConvAEs for the respective variables. Each GCM is passed through the RF, and the median predictions are compiled in Fig. 7A. The mean and mid-50% bounds of current, mid-century, and end-of-century output projections are recorded in Table S4. Ceteris paribus, following model predictions from the tier-1 ScenarioMIP activities, by the 2050s, the median GCM projects grain output to decrease by roughly 107 tons (−0.736% from 2020s) for SSP1, 66 tons (−0.451% from 2020s) for SSP2, 83 tons (−0.568% from 2020s) for SSP3, and 154 tons (−1.055% from 2020s) for SSP5. Compared to RCP-based by-crop studies, this study projects a lower %-decrease in overall grain production up to the 2050s. Specifically, output risk at SSP2 is close to 80% less than that projected for rice by Chen et al. (2020). From the 2050s to the end of the century, grain production would be approximately constant in SSP1, decreasing around 41 tons (107 tons from 2020s, −0.733%) for SSP2, 12 tons (95 tons from 2020s, −0.651%) for SSP3 and 116 tons (270 tons from 2020s, −1.848%) for SSP5. Risk-wise, in the SSP585 scenario, over 50% of models project a lower total grain output in the MLYP provinces than the baseline plotted in Fig. 7A, which is the mean total grain production in MLYP in the first 20 years of the 21st century.

The GCMs can be sorted based on the sum of their point-wise differences based on projected HW frequency (Table S6). Generally, similarities and differences in their AOGCM components impact HW projections of GCMs systemically. The sorted GCMs are then split into two ensembles: the “higher” if their HW frequency projections fall into the top 50%, and the “lower” if their HW frequency projections fall into the bottom 50%. The specific order can be found in Table S6. Figures 7B and 7C split Fig. 7A by GCMs into “lower” and “higher” ensembles according this definition.

Comparing the two ensembles, those projecting higher HW frequencies initially project higher annual grain output. However, as highlighted in ‘HW-CW Relationship’, the more significant decrease in CWs corresponding to higher HW frequencies should be held more accountable for this increase. Moreover, across scenarios, a significantly higher level of output volatility is seen between 2021 and 2060 in the “higher” ensemble, with fluctuations of the 4-year median output of as much as over 150 tons (seen in SSP585 2055-2058 anomaly from 2051-2054). Max output fluctuations in the “lower” ensemble in that same time range is at just over 110 tons, more than 30% less than the “higher” ensemble.

Concerning scenarios, SSP245 grain production projections in the lower ensemble (Fig. 7B) in the second half of the 21st century fluctuated around 14,500 tons and never broke the baseline, but SSP245 in the higher ensemble quickly broke the baseline at 2071-2074, and sustained there, 50 tons below the value the lower ensemble settled around. In the highest radiative forcing scenario, SSP585, the median of MLYP grain output projections by the higher ensemble dives below baseline (and sustains there) as early as the mid-2060s. In comparison, the lower ensemble maintained until the mid-2070s and broke the baseline. Interestingly, starting in the 2060s, with nonstop fluctuations, 4-year median output projection of the SSP370 scenario in the higher ensemble exceeds the output projected of SSP245, and in two 4-year smoothing periods exceeded that of the sustainability scenario, SSP126. However, the median and range of the SSP370 projection do not platform nor slow down in its decrease; instead, regression analysis on the medians shows an increasing slope over time. This indicates that optimistic projections of SSP370 grain output should be subject to further investigation.

Fortunately, as seen across the two ensembles, if sustainable development and land use are achieved, neither of the ensembles projects a substantial decrease in grain output by the end of the century, which, again, is a result explicable by the extra reduction of CWs that neutralized the impact of some heatwaves.

Impact of average temperature change on grain production.

The tas anomaly from the 1981-2010 climatology projected by each GCM are drawn to grain production anomaly in Fig. 8. Adhering to the Central Limit Theorem, anomalies are averaged across each period of interest. Regression analysis excluding outliers shows a clear negative relationship between grain production and tas anomaly across scenarios in the following decades, mid-century, and end-of-century MLYP grain production projections, with regression certainty increasing over time. However, step-wise regression analysis on 2021-2040 points to a roughly 200 tons (1.67%) decrease in grain production with a one degC increase. On 2041-2060, step-wise regression points to an approximately 145 tons (1%) decrease in grain production with a one degC increase. In 2081-2100, grain production decreased approximately 120 tons with a one degC increase in mean tas. Despite the slight, this decrease in response magnitude is noteworthy.

Model performance

According to Table 4, the optimized ConvAE-RF model performs significantly better than the optimized B1-FGLS baseline in all three evaluated metrics. This result aligns with Jeong et al. (2016) and Fradgley et al. (2023), indicating successful model optimization efforts.

The optimized ConvAE-RF model is also more accurate (lower MAPE) and precise (lower EVar) than B2-RF. With the lowest D_γ amongst the three evaluated models, the proposed model is the best fit for the MLYP dataset. It would be safe to conclude that the space-aware RF regression powered by convolutional embeddings performs superior to the common RF regression in the modeling of grain production in the MLYP. This indicates that added spatial complexity can increase the performance of a predictive RF model, which aligns with the observations about RF-based models in adjacent fields, such as Zhao et al. (2022).

However, it should be acknowledged that increased model complexity is making it harder to understand and interpret the results. The black-box nature of neural networks, coupled with the ensemble nature of random forests, complicates the interpretability further. Therefore, discussion of results will put a heavier focus on regional patterns over province-specific results, which may be inaccuarate due to a range of influences (more about such influence in ‘Limitations and Uncertainties’). Additionally, the model performs well in the mid-lower Yangtze plains (MLYP), but its applicability to other regions cannot be guaranteed.

Projected ETEs and their negative influence

In the proposed model, climate-induced decline in grain production in the MLYP is comprised of two general stages: decline and platform. Up to the 2070s, MLYP grain production is projected to decline steadily (with step-wise decline increasing as Scenarios escalate from SSP126 to SSP585) (Fig. 7). The diminishing influence of overall climate change is explained by the mixed influence of ETEs at higher temperature anomalies. This is especially significant in high-radiative-forcing scenarios of SSP370 and SSP585, where CW frequency diminishes significantly as a tradeoff for increased HW. At this stage, production decline from the damaging factors is offset by decreasing CW frequency’s relatively impactful positive influence. Production decrease in the SSP585 scenario lies within the 0%–4% range projected with RCP8.5 by Sun et al. (2023).

Historically, extreme temperatures impact more on grain production than average temperature (Liu et al., 2022). This study supports this conclusion in a timeframe beyond the status quo. From the RF regression, hw_freq’s variable importance on MLYP grain production is 237.4% that of PrimInd, 247.5% that of pr, and 40.5% that of tas; cw_freq’s variable importance on MLYP grain production is 358.7% that of PrimInd, 374.0% that of pr, and 61.3% that of tas (Table 3).

Laboratory studies, such as Rodríguez et al. (2015) and Gao, Tester & Julkowska (2020), indicate that heat stress (appearing in nature as HW events) reduces rates of photosynthetic efficiency and fresh weight observed in grain crops upon exposure in lifetime. Combined with its considerable impact to the performance of the ConvAE-RF model (Table 3), HW frequency can be held fairly accountable for the decrease in grain production. CWs are also capable of hindering rice production (Kim et al., 2021) and negatively impact the quality of crops (Mathivanan, 2021), so significantly decreasing CW frequency over the rest of the 21st century across scenarios (Fig. 4) may be a remedy.

Comparing Figs. 7B and 7C, despite the absence of a solid median output difference, the projected median grain outputs of all scenarios in the higher ensemble are much more volatile than the lower ensemble. Moreover, in the higher ensemble, median curves in SSP126, SSP245, and SSP585 almost mirror the shape of HWs and CWs. This resemblance implies a causality relation, albeit a holistic causality test was not performed due to complexity constraints.

Impact of changing temperatures

RF regression output in Table 3 places mean air temperature (tas) as the 2nd most crucial variable (after total sown area of grain crops, SownArea). In other words, tas significantly influences the model’s performance, and changes in tas are highly associated with changes in MLYP grain production. Number-wise, as in Fig. 8, grain production variability with tas anomaly is most uncertain when tas anomaly is in the range of 2 degC to 3 degC. When tas anomaly is above six degC, none of the models project an increase in grain production from baseline. As in Fig. 7, decreasing production projections under a warming climate corroborates this assessment. However, it should be noted that the influence of temperature change on production in this study seem to be slightly lower than the conclusions of the previous studies, such as Chen et al. (2020) and Holst, Yu & Grün (2013).

Drawing the variability in grain output over time by province provides a plausible explanation (Fig. 9). CW frequency decrease in Jiangxi and Jiangsu and HW & tas increase in northern Anhui (‘Spatial Patterns’) provided more favorable conditions in previously low-temperature crop growing periods, thereby offsetting some negative impacts of climate change in the broader MLYP zone. Province-level output variability (Fig. 9) also warns against significant, negative consequences in high-radiative forcing scenarios, with multiple provinces expected to incur a grain output decline larger 1% in the 20-year periods after 2060. More climate-induced grain production variations is also projected across provinces in higher-radiative-forcing scenarios: a climate-induced provincial grain production fluctuation greater than 1%, whether negative or positive, is never projected in SSP126 (Fig. 9A) or SSP245 (Fig. 9B). For why positive changes may occur (especially Anhui and Jiangxi), please refer to the trade-off between HWs and CWs discussed in ‘Projected ETEs and Their Negative Influence’.

The existing literature does not reach a general agreement on pr’s influence on grain output, regional studies (such as Sharma et al. (2019)) and studies that included the MLYP as a significant part of a specific sector under bigger scales (such as Holst, Yu & Grün (2013)) agreed on the slight negative influence of increasing precipitation in the temperate zone on grain production and a strong negative impact of extreme precipitation (Liu et al., 2022). Yue et al. (2021) claims that precipitation in the future may increase. This provides for a possibility that may participate in neutralizing the negative influence of tas, especially when torrential rains are considered.

Limitations and uncertainties

Due to the high complexity of this model and the diversity of our data sources, there are a number of noteworthy uncertainties in this study.

Model-wise, the combined use of the ConvAE architecture to handle spatial complexity and the RF regression to handle inter-variable relationships has helped this study to perform superior to existing models, as mentioned in ‘Model Performance’. However, the proposed model requires a substantial amount of data and may be less scalable. This could be a limitation to applying the proposed to larger or atypical datasets in different geographical areas and climates.

Data-wise, future climate projections comes from GCMs. While GCM data provides valuable insights into climate patterns, it is subject to inherent uncertainties due to model assumptions, resolution limitations, and the representation of complex physical processes (Thrasher et al., 2022). Consequently, predictions based on GCM data should be interpreted with caution, considering the potential variability and the influence of factors not fully captured in the models. A discussion of how the GCMs in this study exhibit such variations can be found in the Supplementary Section S2. There is also a level of difference in crop type, plant phenotype, and farming techniques across regions. Within the area of study, Anhui and Hubei are the only provinces that predominantly cultivate a single kind of crop. Anhui is a major producer of Paeonia ostii (Peng et al., 2017), while Hubei Oryza sativa (Fan et al., 2018).