Building social resilience in North Korea can mitigate the impacts of climate change on food security

Adaptation based on social resilience is proposed as an effective measure to mitigate hunger and avoid food shocks caused by climate change. But these have not been investigated comprehensively in climate-sensitive regions. North Korea (NK) and its neighbours, South Korea and China, represent three economic levels that provide us with examples for examining climatic risk and quantifying the contribution of social resilience to rice production. Here our data-driven estimates show that climatic factors determined rice biomass changes in NK from 2000 to 2017, and climate extremes triggered reductions in production in 2000 and 2007. If no action is taken, NK will face a higher climatic risk (with continuous high-temperature heatwaves and precipitation extremes) by the 2080s under a high-emissions scenario, when rice biomass and production are expected to decrease by 20.2% and 14.4%, respectively, thereby potentially increasing hunger in NK. Social resilience (agricultural inputs and population development for South Korea; resource use for China) mitigated climate shocks in the past 20 years (2000–2019), even transforming adverse effects into benefits. However, this effect was not significant in NK. Moreover, the contribution of social resilience to food production in the undeveloped region (15.2%) was far below the contribution observed in the developed and developing regions (83.0% and 86.1%, respectively). These findings highlight the importance of social resilience to mitigate the adverse effects of climate change on food security and human hunger and provide necessary quantitative information. The role of social resilience in mitigating hunger related to climate change is explored in North Korea, South Korea and China, regions with similar climatic conditions but varying levels of economic development.

A daptive ability that depends on social-economic resilience is consistently regarded as a crucial factor in coping with climatic risk and safeguarding food security. Social-economic resilience includes not only traditional financial assets and infrastructure 1 but also demographic structure, resource utilization, technology, education, and attitudes and perceptions of risk to change adaptive behaviours 2,3 . However, few studies have integrated social resilience into the climate risk framework to quantify the contribution to regional food production of inherently hard-to-quantify human behaviour and risk perception. Moreover, little is known about the potential of social-economic resilience to mitigate the adverse effects of climate change and extreme events to ensure food security.
North Korea (NK), located in a climate-vulnerable region of eastern Asia (Fig. 1), has been strongly affected by climate change. Many meteorological disasters have induced more severe famine over the past few decades, including typhoons, heavy precipitation events and river floods 4 . According to reports, NK suffered from a freezing disaster in 1993, hail in 1994, severe floods from 1995 to 1996, a typhoon and drought in 1997, and frost in 1998, among other disasters 5 . Even in the twenty-first century, NK's grain production still cannot meet the population's needs, and food deficits still loom large and are even a growing trend 6,7 .
The similar climatic conditions 8 but varying levels of economic development of NK (an undeveloped region according to the World Bank classification) and its neighbours (South Korea (SK), a developed region; and China, a developing region) provide a natural example for investigating the impacts of climate extremes and their link to social-economic resilience 9 . The similar climatic conditions 8 of these areas also rule out uncertainties in social-economic assessments due to differences in climate vulnerability. Comparing economic vulnerability based on social resilience among the three regions with natural adjustment for climate risk is a valuable approach-that is, the ability to adapt to climate risk for food security resulting from climate and economic vulnerability 10,11 . Here we study rice (Oryza sativa L.), as it is one of the most essential foods in NK. Rice composes more than 60% of the total grain production and directly affects food security for NK in terms of planting area and production 12 . Notably, the adverse impacts of climate change on rice systems are increasing 13 .
Obtaining reliable statistics and survey data from NK is difficult due to NK's politics and economics. Therefore, this study attempted to fully use remote sensing and climate data with openly available statistical information to examine and assess climatic risk and food security with the interaction between climate and Building social resilience in North Korea can mitigate the impacts of climate change on food security social-economic vulnerability for NK and its neighbours (Supplementary Table 1 and Supplementary Fig. 1). In addition, the method presented in this study can be used in regions of the world that lack official information to evaluate climatic risk and food security status ( Supplementary Fig. 1). This research intends to answer three interrelated questions: (1) How has climate change (climate extremes) affected rice production in NK in the past? (2) To what extent would projected climate change affect rice production loss in NK in the future? (3) How have human activities (adaptation based on social-economic resilience) exacerbated or ameliorated food deficits in NK and its neighbours? Specifically, we focus on normal and extreme climate changes in NK over a recent past 18-year period (2000-2017), and we attribute rice biomass changes to climatic factors resulting from high-frequency climate extremes and increased vulnerability. The climate projections presented here are based on 27 global climate models (GCMs). The GCMs are derived from the Coupled Model Intercomparison Project Phase 6 (CMIP6) under two future shared socio-economic pathways (SSPs; SSP245 represents SSP2 + RCP4.5, a medium-development pathway; SSP585 represents SSP5 + RCP8.5, a high-development pathway). Moreover, we assess regional climate changes and resulting production losses from a climatic risk perspective. Finally, the effects of social resilience are preliminarily explored on the basis of five factors (population development, resource use, science and education, economic development, and agricultural inputs) to mitigate climate shocks for rice production. Furthermore, the contribution of social resilience to rice production is quantified by contrasting the differences between NK and its neighbours (SK and China). More details on the data used, methods and model robustness checks can be found in the Methods and Supplementary Information.

Results
As the representative of the undeveloped and climate-vulnerable regions, NK depends on rice production to feed approximately 25 million people. Our results show that climatic factors determined rice biomass changes in NK from 2000 to 2017, and high-temperature and precipitation extremes triggered reductions in production in 2000 and 2007. If no action is taken, enhanced climate extremes will cause rice production to decrease by the 2080s, further resulting in hunger in NK. Building adaptation based on social resilience in NK, compared with its neighbours (SK, the developed region; and China, the developing region), can mitigate the impacts of climate change on food security.  Domination of climate over the past 18 years for NK. To attribute changes in rice growth and major climate in NK (Fig. 1), the phenology-and pixel-based paddy rice mapping algorithm was first applied on the Google Earth Engine (GEE) cloud platform to extract rice paddy map referencing used by Dong et al. 14 in which they adequately verified this algorithm and the accuracy of rice maps for ensuring robust application in large-scale fields (Supplementary  Figs. 2 and 3 and Supplementary Table 2). Furthermore, we adopted an ecosystem light use efficiency (eLUE) model to simulate the biomass of the study areas from 2000 to 2017 and calibrated this model. Supplementary Fig. 4 demonstrates that the eLUE model had good robustness and accurately reproduced the LUE of ecosystems monitored by eddy covariance (EC) towers (tenfold cross-validation: coefficient of determination (R 2 ) > 0.75, normalized root mean squared error (nRMSE) < 0.4, P < 0.01). Consequently, we reestablished the distribution of gross primary productivity (GPP) from 2000 to 2017 using this model for NK (Methods and Supplementary Fig. 1). The climatic variables that were screened by variance inflation factors (VIFs) were incorporated in the analysis to determine the climate attribution of rice GPP changes in NK over the past 18 years (2000 to 2017) (Supplementary Table 3). More specifically, for NK, climatic variables explained 80% of the GPP changes observed from flux towers from 2000 to 2017 (the baseline period) (Fig. 2a). In addition, temperature (especially temperature extremes) dominated rice GPP changes (explaining nearly 50%) in NK over the expanded period as seen by the four most important variables (TNn, minimum value of the daily minimum temperature; TR20, count of days when the minimum temperature was >20 °C; FD0, count of days when the minimum temperature was <0 °C; and AAT, average air temperature) being related to temperature (Fig. 2b). These results are robust as indicated in the two validation methods (Fig. 2a and  Supplementary Table 4). Furthermore, in the 18 years from 2000 to 2017, our results show sudden drops in rice production in 2000 and 2007 ( Supplementary Fig. 5). The potential climate shocks induced by climate extremes can cause fluctuations in NK's rice production.
Specifically, extreme heat events and precipitation caused the fluctuations in rice production in 2000 and 2007 (respectively) in NK, as seen by comparing the spatial anomalies of each climatic variable with those in other years. An abnormal increase in extreme heat index was observed in 2000 (TXx, maximum value of the daily maximum temperature in Supplementary Fig. 6), and long-term high temperatures triggered a heatwave that resulted in warming for daytime and overnight periods (TR20 and SU30, count of days when the maximum temperature was >30 °C, in Supplementary Fig. 6). The frequency of abnormal increases in TXx, TR20, SU30 and FD0 in 2000 accounted for 27%, 35%, 37% and 54%, respectively, of the entire region, especially in the western and southern rice-growing areas (Supplementary Figs. 6 and 7). Substantial increases in rain and rainy days were observed in the non-rice part of northeast NK. This did not alleviate the reduction in production caused by the high temperatures and the heatwave ( Supplementary Fig. 6). In 2007, precipitation extremes dominated, decreasing rice production in the west/southwest rice-growing region (TP, total precipitation, and R50, count of days when precipitation was ≥50 mm, in Supplementary Fig. 6). Specifically, the abnormal increases in TP, R50 and R25 (count of days when precipitation was ≥25 mm) accounted for 87%, 72% and 80%, respectively, of the entire region ( Supplementary Fig. 7). Furthermore, the long-term and substantial precipitation produced conditions that made plants highly susceptible to crop root rot and flood damage. These precipitation extremes regulated the surface temperatures, causing the maximum temperatures to decrease and the minimum temperatures to increase (TXx and TNn in Supplementary Fig. 6). Additionally, extreme heat in 2000 and precipitation extremes in 2007 occurred during key phenological stages of rice development-that is, heading-tillering stages that are usually sensitive to climate change and especially sensitive to high temperatures that can cause plant death ( Supplementary Fig. 6). In short, rice production in NK decreased abruptly due to extreme weather events (temporally and spatially), resulting in human hunger. Relative importance (%)   Table 3 for the detailed definitions of TS (total solar radiation), AAT, TP, TNn, TXx, TR20, SU30, FD0, R50, R25, R10 (count of days when precipitation was ≥10 mm) and R1. R 2 CV in b represents the goodness of fit from the fivefold cross-validation (Supplementary Table 4).

Future climate change and production losses in NK.
Here we show the effects of climate change on rice production losses under two climate scenarios (SSP245 and SSP585) based on an ensemble of 27 GCMs (Supplementary Table 5 Fig. 8c,i,j). In general, no matter which climate scenario is considered, the risk of high temperatures and extreme rainfall due to future warming will increase. Consequently, with extremely high temperatures and altered precipitation in the future, rice biomass in NK would decrease by 18.9% under SSP245 and by 20.2% under SSP585 in the 2080s, compared with the baseline period. Production would drop by 13% under SSP245 and by 14.4% under SSP585 in the 2080s (Fig. 3a). In NK, where vulnerability to climate is exceptionally high and the frequency of extreme weather leads to low production, 20.2% biomass losses may be conservatively estimated, and the fragile food system may collapse, resulting in famine. In the future, the negative impact caused by climate change on biomass will extend across NK and show different decreasing trends among the regions. The most serious failure of rice biomass was found in the southwest and on the east coast (Fig. 3b), which are the breadbasket of NK (Fig. 1). Despite sporadic increases in biomass, the overall biomass loss from central NK to western coastal areas hides these inconsistent estimates from different climate models. Rice biomass losses projected by the 27 GCMs will become more consistent in the 2080s under the SSP585 scenario (lower standard deviation), which means that the breadbasket areas will be the most extreme hotspots in terms of decreased rice production ( Supplementary Fig. 9).

Analysis of adaptability in NK, China and SK.
Although rice production is generally subject to natural-environmental change from the standpoint of a climatic risk framework, the adaptive capacity at regional or national levels from social resilience is a greater determinant of rice production losses 15 . Social resilience is driven by population, economics, technology and culture 3 . To quantify the contribution of social resilience to rice production, we collected and used economic statistics from the Food and Agriculture Organization (FAO), the World Bank and an agricultural dataset of remote sensing that involved five factors-population development, resource use, science and education, economic development, and agricultural inputs ( Table 1 and Supplementary Table 6)-that together constituted social resilience. Supplementary Figs. 10-12 indicate the correlations between a single social-economic variable and rice production. NK showed weaker correlations between each of the variables of social resilience than SK and China (ρ < 0.5 and no vital significance). When the climate model incorporated social resilience, more significant explanatory power was shown in SK (P < 0.05) and China (P < 0.01), but this effect was not substantial in NK (Supplementary Table 7). In addition, social resilience controlled more rice production changes for SK and China (P < 0.01), improving the contribution by more than 26% and 100%. Social resilience did not provide additional assistance for NK, which meant that climate dominated rice production changes (Supplementary Table 7). Furthermore, social resilience mitigated climate shocks and even converted adverse effects to benefits in SK and China. Specifically, the notably high temperature and heatwave (TR20, AAT, TXx and SU30) threatened rice production in SK. However, regional nitrogen fertilizer input, rural population and population aged 0-14 years reversed these effects from climate shocks and even promoted higher rice production (Fig. 4). A similar phenomenon was observed in China, where resources for social resilience relieved damage from rainfall extremes. Specifically, the changes in access to electricity mitigated the negative influence caused by rainfall extremes (R25) (Fig. 4). This significant moderating effect was not observed in NK, where all interactions between climate shocks and social resilience did not result in increased production (Fig. 4).
We conducted a more comprehensive analysis on the basis of random forest (RF) regression of economics (RF e ) to assess the differences in the contributions of social resilience to rice production among NK, SK and China. Of the 12 indexes of social resilience, higher education, rural population and population ages 0-14 dominated rice production variations in NK (P < 0.05). Patent applications, population ages 0-14, gross domestic product (GDP) per capita and energy use contributed to rice production variations in SK (P < 0.01). The variation in rice production in China was mainly determined by population ages 0-14, rural population, net official development assistance (ODA) received per capita and GDP per capita (P < 0.05) (Fig. 5a). Economic development and population structure in both developed and developing regions played a more critical role than in undeveloped regions, and science and education in the developed area had a greater influence on rice production than other factors (Fig. 5a). The RF e analysis further illustrated the essential contribution of social resilience to rice production in SK and China (explaining 83.0% and 86.1%, respectively, of interannual rice production variation, P < 0.05), which was much higher than in NK (15.2%, P < 0.05) (Methods and Fig. 5a). Rice production in developed and developing regions is therefore controlled by social-economic factors (Fig. 5a). In addition, the robust results obtained in this study demonstrate the great potential of social resilience to increase crop production and resist the harmful effects of climate shocks and weather extremes on food security.
Rice production in NK exhibited a threshold-like response to all four important variables influencing rice production, including school enrolment, tertiary (Fig. 5b). Rice production declined when school enrolment was higher than −0.6 units, presumably because limited capital was invested in science education and reduced economic investment and labour from family sources in agriculture. The nonlinear response of rice production in China to population development provided coherent evidence (Fig. 5b). Specifically, rice production decreased with the increasing agricultural population and population ages 0-14 in the undeveloped and developing regions (NK and China) (Fig. 5b). The pressures from dietary needs caused by population growth and the uncoordinated structure of the population constrain economic development and production increases in undeveloped regions 16 . For instance, the rural areas with relatively higher mechanization contribute to rice production and require lower agricultural populations 17 . In contrast, the response of rice production to population ages 0-14 was the opposite in the developed region (SK), which might be explained by differences in the structure of the population among districts. Developed regions need to increase the proportion of adolescents in the population to adjust for serious ageing 18 so that an abundant agricultural labour force is available to increase rice production. Energy use produced an increase in rice production for the developing region (China). Yet, the opposite result was observed for the developed region (SK), owing to the capacity to import resources because of sufficient capital 19 . Social resilience probably results in various impacts on rice production among NK, China and SK.
The nonlinear response of rice production to agricultural inputs in the three regions is shown in Supplementary Fig. 13. The response curves for nitrogen, phosphorus and irrigation showed increasing trends in SK and China. However, in NK, the production responses to the three agricultural practice inputs were expressed as humps or concave curves ( Supplementary Fig. 13). On the basis of the RF model and its out-of-bag error, agricultural inputs (nitrogen, phosphorus and irrigation) explained −4.8%, 51% and 77% of rice production changes for NK, SK and China, respectively. The explanatory degree for NK was negative, and the three agricultural practices cannot support the increase in rice production in the current situation.

Discussion
Our results demonstrate that climate extremes reduced rice production and that climate warming would contribute to famine in NK. The high-temperature index (SU30) would increase 97.6% and 221.94% under SSP245 and SSP585, respectively, by the 2080s (Supplementary Fig. 8). The number of high-temperature days per year would increase by nearly one month by the 2080s under SSP245, and the number would increase by two months under SSP585. This finding is also supported by results documented in previous literature. Kawasaki and Uchida 20 indicated that the abnormal temperature caused by global climate change had a severe negative impact on agricultural production and quality. And climate extremes will increase in frequency and intensity with rising temperature in the future 21 . Our results also suggest that projected GPP and production may decrease by 20.2% and 14.4% in NK in the 2080s under the SSP585 scenario. However, this result may be underestimated and could be even more severe for crops, because we used the average GPP and harvest index for entire regions to calculate relative changes in biomass. Climate extremes have been known to have critical impacts on the resilience of the food supply chain 22 . Temperature extremes observed in the past have contributed to increased yield variability, and extreme temperature events will probably increase in the future 23 .
However, abnormal weather and climate extremes do not always have negative consequences. Some studies relating to differing weather conditions and intensity and target crops have reported different results in various areas. For instance, Zhao et al. 24 indicated that rice yields declined significantly with higher temperature on the basis of simulations of field warming experiments, statistical models and gridded crop models. These yield reductions exceeded assessments by the International Food Policy Research Institute 24 . Yet, in India, rice yield in the most agro-ecological zones would benefit from climate change predicted by climate scenarios and GCMs 25 . Furthermore, Lesk et al. 26   rainfall (>50 mm hr −1 ) caused severe damage to crop yield but that crop growth benefited from heavy rain of 20 mm hr −1 . Due to the projected increases in temperature and rainfall intensification in the future, the impact of climate extremes on crops remains somewhat uncertain. It is difficult to identify major changes in climatic sensitivity solely on the basis of production data over time because of the infrequency of extreme weather 27 . Inductive analysis of the changes in exposure and crop sensitivity to climate extremes is therefore a prerequisite for food security assessment. Differences in the quantitative attribution in rice production between NK and its neighbours (SK and China) reflected the importance of adaptation based on social resilience to mitigate the adverse effects of extreme climate. Food security may benefit most from changes in adaptive capacity under future climate change, such as agricultural practices, economic development 28 , resource use and social cognition 29 . Irrigation 30 , fertilization, conservation tillage 31 and crop breeding are being given much attention, with more focus on increasing yields. For example, Challinor et al. 32 found that anthropogenic adaptability increases the average yield of crops by 7-15%. However, suppose the focus is only on increasing crop yields. In that case, farmers will still suffer severe production losses due to their low-risk perception concerning the effects of extreme weather disturbances. Learning reflects the ability to produce, absorb and transform new information about climate risk, adaptation and coping with uncertainty. This ability to learn and apply further scientific knowledge is a mitigation mechanism applicable to climate change 33 . Capital investment that depends on economic development is a more direct approach. These investments include building early warning systems and climate insurance, resource and energy use, and international trade, as well as reducing poverty. Insurance is a tool to mitigate climatic risk and restore livelihoods, especially in response to climate extremes 34 . Still, if the insurance structure is not correct, it has an inhibiting effect on risk reduction 35 . For regions that cannot compensate for losses through trade, these years of low productivity can still be devastating. Undeveloped countries are more vulnerable and less resilient to climate change.
When the poor are struck, they have less support from friends, family and the financial system. Policies meant to reduce poverty under similar climatic risk conditions can also reduce the adverse impacts of climate change 36,37 . The interaction of social-economic factors from many aspects is the key to decreasing economic vulnerability and increasing social resilience to mitigate climate vulnerability. In this study, social resilience was shown to be enormously important for reducing hunger by contrasting the situation in three regions. Moreover, the topography may limit future rice production potential in NK due to the poor-quality land that is not suitable for planting rice (Supplementary Section 1.1 and Supplementary Fig. 18). Future planning for food security needs to consider climate change and social-economic interactions and development.
This study provides critical insights into the contribution of social resilience to food security. The regions of food insecurity in the world are also the regions with limitations to necessary data due to conflicts, war and climate extremes. A comprehensive evaluation of the food security status in these regions is not only an imperative requirement for sustainable development but also a necessary way to decrease hunger. This study integrated meteorological reanalysis products, openly available statistics and high-resolution imagery to replace unavailable data and ensure robustness. These analyses of food production and food security can therefore be extended to other parts of the globe with ground-truth data limitations. The platform employed in this study, coupled with distinctive regional characteristics, can analyse the importance of social resilience to food security. It has the potential to produce comprehensive assessments of food security status in other regions of the world, such as Africa, Latin America and Southeast Asia. Predictably, national or regional food security is more attributable to social-economic features. Given that food loss occurs from different causes and risks beyond the climate risk framework, global food production is subject to the constraints imposed by exposure and adaptation. For example, food insecurity in most areas of Africa is attributable to an increase in violent conflict since 2014. However, production losses are also attributable to locusts and drought. New insights regarding food security can be obtained for specified regions by incorporating different factors into exposure. This will enable us to further explore the ability of social resilience to mitigate factors that destabilize food security. With the future development of multivariate data (especially remote sensing data), study regions lacking reliable data will have new opportunities for analysis and evaluation.
Our findings are subject to some limitations and should be considered with caution. First, the lack of sufficient ground-truth data for vulnerable regions was the biggest challenge, particularly socialeconomic data (only 20 years). Additionally, we do not resolve other potential sources of uncertainty. For example, the lack of short-term fertilization data and additional management data is a potential source of uncertainty not considered in this study 38,39 . Second, soil properties may also introduce uncertainty in the results 40 . The physical and chemical properties of soil were different in each region, yet we do not include these in our model. Folberth et al. 41 found that assessing the impact of climate change on yield depended on soil type because soil characteristics and moisture buffer or amplify climatic impacts 27 . This study also did not consider the effects of rice genotype due to a lack of cultivar data, and therefore there may be uncertainty regarding regional production differences 42 . For GPP loss estimation, we did not introduce carbon dioxide (CO 2 ) concentration in our model due to the controversy regarding the ability of increasing CO 2 concentration to decrease crop water demand and improve yield while reducing the nutritional content of grain 43 . Not introducing CO 2 concentration effects may have affected reliable and robust climate and economic vulnerability estimates. Finally, the rice maps and downscaling methods were also sources of uncertainty. Specifically, uncertainties regarding rice maps may have arisen from remote sensing data, poor weather, data quality     Table 1 and Supplementary Table 6 (Table 1). *P < 0.1; **P < 0.05; ***P < 0.01. R 2 indicates the explanation degree of the nonlinear model for contribution of social resilience to rice production in the three regions. MSE, mean squared error. b, The nonlinear response of rice production to standardized values of the four most important variables of social resilience. The black lines are smoothed representations of the response using best-fit polynomial equations, with fitted values (model predictions) for the calibration data. The trends of the lines, rather than the actual values, describe the nature of the dependence of rice production on social resilience. The shaded bands denote the 95% confidence intervals. The green and orange dashed lines in the partial dependence plots represent rice production from the baseline and SSP585 scenarios, respectively. and the number of images observed 14,44 . This study generated future daily meteorological data with statistical downscaling, but different downscaling methods (such as the change factor method and dynamical downscaling) produced different results 45 . In the future, NK's single-crop production is more likely to lead to food insecurity (rice production accounts for more than 60% of the total food production), leading to hunger. Switching agricultural structures and increasing crop diversity may therefore be a critical path to ensuring food security. The introduction of drought-tolerant and heat-tolerant crops is expected to contribute to food security in vulnerable regions such as NK. Especially at the national level, increasing crop diversity by planting multiple crop species can effectively reduce the losses in national food production, thus ensuring the overall food supply 46 . As one of the crucial components of food system resilience, crop diversity can drive trade diversity, dietary diversity and ecosystem services 47 . Climate change is expected to increase the frequency and intensity of climate extremes, probably reducing global food production and famines. An accurate assessment of food insecurity must be valued in the deprived areas of the world because it is a vital link in the world food system and is an essential component of the United Nations' Millennium Development Goals. NK is representative of the world's deprived areas, and it and its neighbours offer an excellent example for studying food insecurity in undeveloped regions. Our robust results suggest that social resilience provides a reliable path for understanding and mitigating the detrimental effects of extreme weather shocks in the future. Using this method, we can clarify future food risks and provide quantifiable pathways and goals for efforts. That will contribute to improved national risk awareness, make up for weaknesses in food programmes and guide adjustments to food strategies optimizing social-economic policies.

Methods
Workflow. The workflow presented in this study was designed to assess the contribution of climate change (climate extremes) and social resilience to regional food security. We can use flexible multisource data to extend the study to the rest of the globe, where data availability is restricted (Supplementary Fig. 1). Regional biomass was first simulated on the basis of an eLUE model using observations of flux towers and a remote sensing index (with a 500 m spatial resolution and an eight-day (8d) time step) from NK's neighbours (Liaoning and Jilin provinces of China (CHN_1_2)). We also conducted cross-validation for predicted biomass. Second, we calculated 12 normal and extreme weather variables from ERA-5 reanalysis (with the original 0.1° spatial resolution and a daily step) for NK, SK and CHN_1_2 during the rice growing seasons in 2000 to 2017. These weather variables were used as predictor variables to build a RF model. Furthermore, annual rice distribution was extracted by phenology-and pixel-based paddy rice mapping and satellite products from the Moderate Resolution Imaging Spectroradiometer (MODIS) to mask non-rice regions.
Then, for climate shocks and future climatic risk, two steps were employed. In Step 1, gridded climate variables and biomass in the rice regions from 2000 to 2017 were used as the baseline to build the RF model (RF b ). We also validated the stability of the model and calculated the importance of the variables. In Step 2, NK's climatic risk and biomass losses were projected on the basis of the model determined in Step 1 with future climate data that were used for statistical downscaling with daily climate data from ERA-5 (1979 to 2017) and 27 GCMs under SSP245 and SSP585 to assess NK's food security status by the 2040s and the 2080s. Finally, this study included social-economic variables (derived from FAO statistics) and climatic variables (based on daily climatic products of ERA-5) at the national level from 2000 to 2019. The potential for social resilience to mitigate climate shocks in NK, SK and China was explored using a moderation model. We also focused on nonlinear responses of critical variables affecting rice production using an RF e model and conducted time-series split cross-validation to further clarify the contribution of social resilience to rice production and mitigating climate shocks.
Data sources. We used daily reanalysis data to analyse climate variables over the years of the study. The daily reanalysis data were obtained from the European Center for Medium-Range Weather Forecasts' (ECMWF) gridded dataset at 0.1° resolution. This dataset included daily 2-m temperature (24-hour maximum, minimum and mean temperatures), precipitation and solar radiation from 1979 to 2018. See Supplementary Table 1 for more details. We selected the ECMWF's ERA-5 dataset for two reasons: (1) ground data are not readily available, and (2) by using the global dataset, we can apply our methods to other areas with limited data (Supplementary Table 1).
Our primary analysis calculated three average and nine extreme climate indexes for the 1979-2017 and 2021-2100 periods (Supplementary Table 3). Climatic variables for 2000-2017 were used to attribute extreme weather, modelling regression and prediction for NK. We projected future climatic risk and production losses under different climate scenarios using variables from 2021-2100. We also conducted statistical downscaling of the data from 1979 to 2017 to project future climate change using ERA-5.
We employed the MODIS gridded mosaic for remote sensing information, including 8d surface reflectance, 8d leaf area index and 8d GPP accessed from GEE for 2000-2017 (Supplementary Table 1). We used MODIS surface reflectance to calculate vegetation indexes-that is, normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), land-surface water index (LSWI) and the near-infrared reflectance of vegetation (NIR v ). The specific formulas are: where ρ Red , ρ NIR , ρ SWIR and ρ Blue are the surface reflectance values of the red band, the near-infrared band, the shortwave-infrared band and the blue band in the MODIS imagery; L is the canopy background adjustment that addresses nonlinear, differential near-infrared and red radiant transfer through a canopy; C 1 and C 2 are the coefficients of the aerosol resistance term, which uses the blue band to correct for aerosol influences in the red band 48 . And G represents the gain factor. The parameters in the EVI formula are, L = 1, C 1 = 6, C 2 = 7.5, and G = 2.5. For geographic information, we also used the digital elevation model with 90-m spatial resolution from SRTM Digital Elevation Data Version 4 in the GEE platform to calculate the slope for every grid cell. All datasets supporting the results of this paper are freely available in Supplementary Table 1.
We accessed the observations of daily net ecosystem CO 2 exchange and ecosystem respiration from two EC towers in the Chinese FLUX Observation and Research Network from 2003-2010 to calibrate and simulate the biomass of the study areas without meteorological inputs (except solar radiation). The EC towers were located close to NK in Yucheng, Shandong province of China (116° 34′ 12.72″ E, 36° 49′ 44.4″ N) and at Changbai Mountain, Jilin province of China (128° 5′ 45″ E, 42° 24′ 9″ N), with farmland and forest ecosystems, respectively. We used the daily net ecosystem exchange and heterotrophic respiration from the flux towers to calculate the gross ecosystem CO 2 exchange. We refer to gross ecosystem CO 2 exchange as gross ecosystem primary productivity (Supplementary Table 1 44 . The reason for choosing these sites was based on the similar climate zones in the study areas (Fig. 1) and the lack of actual data from NK. All data were strictly examined to meet the standards for further analysis, including cross-validation and comparison with existing data.
The statistical data were obtained from FAO (rice production and population), the United Nations Statistics Division (GDP and imports/exports of goods and services) and the World Bank (population ages 0-14, population ages 15-64, rural population, energy use, access to electricity, school enrolment, patent applications and net ODA received per capita) ( Table 1 and Supplementary  Table 6). All data are based on the period 2000-2019 at the country level. Nitrogen and phosphorus fertilizer applications for agriculture production were provided by Lu and Tian 49 , who developed the global gridded data at 0.5° × 0.5° resolution from 1961 to 2013 and published for free access. The statistical data were used to analyse the social-economic attributions. FAO's data quality is generally divided into ten categories: 'Unofficial figure' , 'Symbol for indigenous or liveweight meat' , 'Official data' , ' Aggregate (may include official; semi-official; estimated or calculated data)' , 'Calculated' , 'FAO estimate' , 'Calculated data' , 'FAO data based on imputation methodology' , 'Data not available' and 'Trend' . More than 95% of the data used in this study were from the 'Official data' and ' Aggregate' categories to ensure high quality.
For further analysis of non-climate attributions, we considered the effect of irrigation on rice growth. We used the water consumption coefficient for rice paddies to replace irrigation because of the lack of irrigation data over the study areas and the large uncertainty in irrigation timing, amounts and methods. The irrigation period for these maps was from 2001 to 2017, and we calculated the ratio of evapotranspiration (from MOD16A2, Supplementary Table 1) to precipitation as the water consumption coefficient. Furthermore, we averaged the water consumption coefficient for rice paddies in the study areas.
Rice paddy mapping and estimating biomass. Annual rice maps were used to mask non-rice regions to improve the accuracy of the regression model. The annual gridded rice paddy maps (from 2000 to 2017) were produced on the GEE cloud platform on the basis of vegetation indexes (equations (1)-(3)) from MODIS reflectance bands and the calculated rice phenological calendar (Supplementary  Table 8 and Supplementary Figs. 2 and 3). See Dong et al. 14 for clear steps and complete verification. We used very high-resolution samples from Google Earth images to validate the 2015 rice map and obtain high-accuracy rice distributions over NK, SK and CHN_1_2 (Supplementary Table 2 and Supplementary Fig. 3).
Existing biomass products are usually calculated using meteorological data. However, a misleading result can be obtained when analysing which climate variables dominate biomass changes using these products-that is, the climate variables that dominate biomass changes are highly correlated with the climate factors used for the calculation. We therefore considered using a vegetation index from MODIS (with a 500 m spatial resolution and an 8d step) with daily EC tower observations from NK's neighbours (CHN_1_2) to estimate 8d GPP without extra climatic factors (such as temperature and precipitation) for NK. Traditional LUE models need to assess the fraction of absorbed photosynthetically active radiation (fAPAR) and actual LUE (ε) separately when calculating GPP. Furthermore, using VIs × PAR TOC to evaluate GPP means that vegetation indexes (VIs) can be more clearly used as a measure of eLUE (described as LUE based on photosynthetically active radiation (PAR)) 50 . eLUE is defined as the ratio between GPP and PAR at top-of-canopy (PAR TOC ): We averaged daily GPP (GPP EC , g C m −2 d −1 ) and daily PAR TOC (MJ m −2 d −1 ) from two EC towers as 8d values to calculate eLUE TOC (g C MJ −1 ). f(VIs) was from the regression of eLUE TOC to VIs. Once eLUE TOC was estimated, the eLUE relationship with PAR TOC (equation (5)) was rearranged to predict GPP: f(VIs) is the coupling relationship between eLUE and VIs. Past research has focused on the linear model 50,51 . However, multiple variables have nonlinear relationships in the natural environment, and large-scale areas contain multiple complex ecosystems. The performance of the linear model can no longer meet the practical application in large-scale regions or be extended to other districts. We converted f(VIs) into a nonlinear model (that is, an RF model) and incorporated a variety of VIs into the model (NDVI, EVI, leaf area index and NIRv). RF models have good performance and fewer parameters than other nonlinear models. In recent years, RF models have been widely used in different regions to solve natural science problems at the global scale 52 . For a more detailed explanation of RF models, see Breiman 53 .
To establish the relationship between eLUE and VIs (calibrating the eLUE model) and provide independent verification, we randomized EC tower data (397 samples) and then divided the data into two subsets for calibration and validation datasets. To evaluate nonlinear model performance, we used a stratified tenfold cross-validation. We used R 2 CV and nRMSE CV to evaluate the results of the cross-validation, where the CV subscript represents the data obtained from the cross-validation datasets 45,52 . We used the Cal subscript to represent the data obtained from the calibration datasets. These parameters were calculated as: where y(i) and x(i) are the simulated and observed values, respectively; y m and x m represent the mean value, respectively, from simulated and observed series; n is the number of samples; pe i and p i are the observed and simulated values, respectively; and pei is the mean value of the observations.
Attributing predominant climate variables in NK. We downscaled 12 climatic variables (originally 0.1°) to 500 m resolution and calculated biomass with 500 m resolution for regional GPP from 2000 to 2017 over NK (see Supplementary   Table 2 for details on the climatic variables). The climatic variables and biomass were then incorporated into an RF model to determine the predominant climate variables. Specifically, all gridded climatic variables with a resolution of 0.1° pixels were resampled to 500 m resolution using the bilinear interpolation method to match biomass. This study's complex models and calculations utilized a large amount of gridded data, resulting in a heavy computation burden and longer times. Hence, we removed the climatic variables whose VIF was greater than 10 to check the multi-collinearity and reduce the computational load and interference in the attribution analysis. The VIF formula was: where R 2 represents the coefficient of determination between the ith independent variable and other independent variables. Therefore, the variables for solar radiation and rain were excluded from the RF b model.
Compared with linear regression, the nonlinear model explained the nonlinear responses of the climate variables and unravelled the influence of related variables. The process used rice maps to mask non-rice areas. The RF model generated different regression trees using random multiple training sets and features, and each regression tree was sampled independently and distributed identically 53 . Each regression tree produced different results through branching, and the prediction from the RF regression model was the average of all trees 53 . The unbiased estimate of RF performance was from out-of-bag (OOB) error, and this result was similar to k-fold cross-validation. For more details about the RF model, please see Breiman 53 and Shi et al. 52 . Specifically, RF b models (Supplementary Table 9) were used to determine the fit between biomass and climatic variables with the parameters 'mtry, the square root of the variables; ntree, 500' in the R program (version 4.0.2; available at https://www.r-project.org) randomForest package (version 4.6-14; available at https://cran.r-project.org/web/packages/randomForest). We further used function importance to compute the variable contribution. The implication of 'importance' is that the permutation of a variable changes the degree of model accuracy-that is, the error rate calculated from samples that did not participate in the training, usually called OOB error. In this study, an increase in MSE percentage was adopted to measure the importance of a variable. Specifically, a variable is permuted, but other variables of OOB samples are kept constant. The RF model is then rerun to obtain new permuted results from OOB. The importance of variables is the average difference between the permuted OOB samples and the original OOB 54 . Mathematically, the importance of explanatory variables is defined as: where y 2,i is calculated from the original OOB samples of RF, y 1,i is calculated from the permuted OOB samples and N tree is the number of trees of the RF model. Commonly, larger VI indicates that the explanatory variables are more important. For NK, SK and CHN_1_2, the percentage of variance explained was calculated by OOB of the RF model as a goodness of fit to assess the response of the climatic variables to biomass (equation (13)) 55 . We validated the stability and accuracy of models on the basis of fivefold cross-validation (Supplementary Section 1.2): where δ 2 y was computed with n as the divisor (rather than n − 1). MSE OOB is the mean of squared residuals, calculated as: where ŷ OOB i is the average of the OOB predictions for the ith observation. We calculated the bias between a specific year and other years during 2000 to 2017 in NK. The ratio of the difference between climatic variables from a specific year and the mean value from other years to the standard deviation of the climatic variables of other years was used to determine the climate anomaly for a specific year. If the anomaly was >1 or <−1, we considered that the values were clearly greater than or less than the others 56 . Specifically, we used the threshold of one standard deviation by assuming that the variations of climatic factors under normal conditions were generally located in the range of one standard deviation about the multi-year mean 54 .
Projections of future climate and biomass losses. We used daily maximum and minimum temperatures, precipitation, and solar radiation from 1979 to 2017 with statistical downscaling and 27 GCMs to calculate climatic variables for different climate scenarios in the future (2021 to 2100), and further calculated three average and nine extreme variables (Supplementary Table 3). The RF b (Supplementary Table 9) was then used to project annual rice losses under future climate change conditions for spatial pixels over NK from 2000 to 2017. Specifically, we used the statistical downscaling model NWAI-WG 57 to downscale monthly gridded data from GCMs to daily climate data for 1,299 grid points from meteorological reanalysis over NK. Statistical downscaling consisted of three major components: spatial downscaling, bias correction and temporal downscaling. Spatial downscaling used inverse distance-weighted interpolation based on the centre of the nearest four grid points in GCMs to improve accuracy 57 , and we then applied bias correction to generate bias-corrected monthly data using a relationship between the observations and GCM data for a historical training period (in this case, 1979-2018). Finally, the daily time series for maximum and minimum temperatures, precipitation, and solar radiation for each pixel were downscaled from the bias-corrected monthly GCM projections using a modified version of the stochastic weather generator 58 . For a more detailed description of statistical downscaling, please see Liu and Zuo 57 .
We downscaled 27 GCMs of CMIP6 (Supplementary Table 5) using the statistical downscaling method as two future climate scenarios (SSP245 and SSP585). Furthermore, the daily average and extreme climate variables were calculated for two climate scenarios from 2021 to 2100 according to the climatic index described above in Supplementary Table 3.
We ran the RF b model on the basis of gridded climate variables and biomass over NK to project future annual biomass and production losses. Before this, multi-collinearity analysis was conducted on the climate variables (equation (11)). More details about the RF model appear in the previous sections. Two validations were used to examine the robustness of the regression model. First, 75% of the random climate and biomass data from 2000 to 2017 were used as training data, and 25% of the data were used as validation data to test the model performance. The R 2 values were used to evaluate the results of validation and calibration (equation (7)). Second, fivefold cross-validation was conducted (Supplementary Section 1.2 and Supplementary Table 4).
We finally assessed climate change and production loss in the future for the 2040s (2021-2060) and the 2080s (2061-2100). Specifically, for climate change in the future, the relative changes of temperature (AAT, TXx and TNn) were calculated by subtracting the means of the historical period (1979-2017) from the future temperature. For biomass, TS, TP, TR20, SU30, FD0, R50, R25, R10 and R1 (see Supplementary Table 3 for the detailed definitions), the relative changes were derived from the ratio of future means to historical means. For production losses in the future, we obtained the inferred mean value of the conversion coefficient (noted as α) from FAO statistical production (production F ) and estimated regional biomass (that is, α as the ratio of biomass to production). The formulas are: where P sta is statistical production from FAO, ∑ B t1 represents the sum of each gridded biomass (GPP projected by the RF b model), α i represents the conversion coefficient for every year, a mean is the average conversion coefficient calculated by α i for every year (2000 to 2017), n is the number of years and A is the pixel area (500 m × 500 m).
Future production was projected by using the mean conversion coefficient and predicted GPP by the RF b model: where P pre is future production and ∑ B t2 represents the sum of each predicted GPP. Finally, we used future and historical GPP and production to calculate relative changes as the losses.
The contribution of social resilience to rice production. For NK, SK and China, we examined the contribution of social resilience to mitigating climate shocks and rice production based on national-level climate variables (Supplementary Table 3) and social-economic variables (Table 1 and Supplementary Table 6) from 2000 to 2019. In addition, the nonlinear responses of the four crucial variables to rice production F was explored by the partial dependence for NK, SK and China to describe how social resilience affected rice production changes.
The selection of these variables was based on evidence from previous studies demonstrating that social resilience depends on these economic and social characteristics. A vital understanding of resilience is that the system may transition to another state once a threshold is exceeded, and it is difficult to return to the original state. Thus, resilience is generally regarded as the ability to persist or absorb change while maintaining the same structure and function, and it is described as the magnitude of system change that can be absorbed or that will result in disruption. Social resilience is based on the interdependent relationship between human activities and ecosystems. People are generally more adaptive to social change when they have access to various financial, technological and service resources and to assets. Asset-based societies are usually more adaptive to environmental changes than poor societies, and the rich are more resilient than the poor. Education in social resilience is the key to explaining increasing productivity and income, and a large part of the demographic dividend is derived from the education dividend. In addition, changing population structure is a direct socio-economic challenge for predicting mitigation and adaptation to inevitable climate change, in which ageing and labour force changes are identified as fundamental socio-economic issues. Energy is a fundamental factor in socialeconomic development and its ability to eradicate poverty. The preamble of the United Nations' Millennium Development Goals recommends "universal access to affordable, reliable and sustainable energy" and recognizes that "social and economic development depends on the sustainable management of the earth's natural resources" (https://sdgs.un.org/goals/goal7). This study further considered necessary management practices for crop production to resist external stress. The relevant variable references are displayed in Supplementary Table 6. The social resilience data were characterized as one of two types to fully reveal the real social situation: soft-adaptive and hard-adaptive 59 (Table 1).
We interpolated missing statistics from FAO and the World Bank (national level) from 2000 to 2019 on the basis of linear regression due to discontinuous economic data (four social-economic variables: energy use, school enrolment, access to electricity and patent applications). The range and the filed and interpolated results of the missing data are shown in Supplementary Fig. 14. Specifically, discontinuous economic data for NK were interpolated using regression models to estimate the missing values (energy use, school enrolment, access to electricity and patent applications) (Supplementary Fig. 14; the R 2 values are 0.70, 0.99, 0.99 and 0.92, respectively; the P values are all <0.01).
A moderation model was used to explore the effects of social resilience on mitigating climate shocks 60 . First, we averaged all gridded climate variables (Supplementary Table 3) for NK, SK and China (here using the climatic data in CHN_1_2 as the replacement) to obtain the region-scale value and match the socio-economic variables. Second, all variables were normalized to better compare the positive and negative effects of rice production changes. Furthermore, we screened all adverse climate shocks and positive social resilience to rice production changes on the basis of the linear regression model. The screened negative climate shocks (as independent variables) and positive social resilience (as moderators) were used for pairwise combinations to build a potential production moderation model. Finally, we introduced the interaction between independent variables and moderators in the potential models to explore the mitigation effects of social resilience. Specifically, if the relationship (the direction and size of the slope of the regression) between two variables of interest (the dependent, Y, and independent variables, X) depends on a third moderating variable (the moderator, Z), moderation is said to occur. Here we considered the simple moderating modelthat is, the following relationship-for hypothesis testing: where Y, X, Z and XZ represent the dependent variables (rice production), the independent variables (negative climate shocks), the moderator (positive social resilience) and the moderating terms, respectively; a, b and c are the coefficients of each regression term; and ε is error. If coefficient c is significant (P < 0.1), then Z represents a significant moderating effect. Analysis of variance was used to examine the significance of different combinations of variables. Given that the effects of social resilience on rice production F are often nonlinear, RF e was expected to perform well in assessing the nonlinear relationship. This study modelled the relationship between social-economic variables and rice production F on the basis of the RF e model (Supplementary Table 9) and tested significance for a single variable and the full model using the A3 package (reference, version 1.0.0; available at https://cran.r-project.org/web/packages/A3) in R. Details about the RF regression model and its parameters were provided in the previous section. The parameters of nonlinear models (ntree and mtry) were set as 500 and the number of the square root of the variables, respectively. The increase in MSE percentage and variance explained were further used as indexes to evaluate variable importance and goodness of fit measures for rice production F (see equations (8), (13) and (14) and the previous section). The partial dependence (that is, marginal effect) was constructed to assess the nonlinear response between rice production F and each of the first four variables. This process was accomplished using the partialPlot function of the randomForest package in R. The threefold time-series split cross-validation was applied to examine the robustness of the RF e model for NK, SK and China (Supplementary Table 10). We divided the training and validation sets along with time series compared with the original cross-validation. The data from 2000 to 2012 were a test of the first fold, and the data from 2013 to 2014 were the validation of the first fold; the data were further expanded at the two-year step. https://srtm.csi.cgiar.org/. The 27 downscaled GCMs of CMIP6 were provided by D.L.L. and H. Zuo 57 , who downscaled them on the basis of the original CMIP6 at https://esgf-node.llnl.gov/projects/cmip6/. The statistical data are freely available from the following sources: data on rice production, rice imports and exports, fertilizer application, and population from FAO are at http://www.fao.org/faostat/ en/; data on GDP from the United Nations Statistics Division are at https://unstats. un.org/unsd/snaama/Basic; and data on population ages 0-14, population ages 15-64, rural population, energy use, access to electricity, school enrolment, patent applications and net ODA received per capita from the World Bank are at https://data.worldbank.org/. Source data are provided with this paper.

Code availability
The first and corresponding authors are prepared to respond to reasonable requests for code.