Mapping the future afforestation distribution of China constrained by a national afforestation plan and climate change

Afforestation has been considered a critical nature-based solution to mitigate global warming. China has announced an ambitious afforestation plan covering an area of 73.78×10⁴ km² for the period 2020–2050. However, it is unclear which areas will be suitable for afforestation under future climate change. Here, we carried out a finer-resolution (25×25 km) dynamical downscaling of climate change for China using the Weather Research and Forecast (WRF) model nested with the bias-corrected MPI-ESM1-2-HR model. Then, using the Holdridge life zone model forced by the WRF model output, we mapped the climatological suitability for forests in China. The results showed that the potential forestation domain (PFD) at present (1995–2014) approximated 500.75×10⁴ km², and it would increase by about 3.49 % to 518.25×10⁴ km² in the period 2041–2060 under the Shared Socioeconomic Pathway (SSP) scenario (SSP2-4.5). Considering the expansion of the future PFD due to climate change, the afforestation area for each province was allocated to grid cells following the climatological suitability for forests. The new afforestation grid cells would be located around and to the east of the Hu Line (a geographical division stretching from Heihe to Tengchong). Due to afforestation, the land cover would be modified. The conversion from grasslands to deciduous broadleaf forests in northern China took up the most area, accounting for 40 % of the new afforestation area. The grid-cell-resolved afforestation dataset was consistent with the provincial afforestation plan and the future climatological forest suitability. The dataset would be valuable for investigating the impacts of future afforestation on various aspects, including the carbon budget, ecosystem services, water resources, and surface hydroclimate regime.

1 Introduction

Afforestation has been considered a reasonable nature-based solution for global warming (Rohatyn et al., 2022; Yu et al., 2022). Afforestation could increase carbon stocks in terrestrial ecosystems by absorbing atmospheric carbon dioxide through its biogeochemical effect (Jayakrishnan and Bala, 2023; Zhu et al., 2019; Gundersen et al., 2021). Moreover, afforestation changes the surface energy and mass budgets as well as the water cycle by modifying the surface albedo and roughness along with the partitioning between sensible and latent heat fluxes (Bonan, 2008; Breil et al., 2021; Wang et al., 2023). Specifically, afforestation causes warming effects through decreasing albedo and cooling effects through intensifying evapotranspiration, which can partly offset or amplify the cooling effects due to carbon uptake from the atmosphere (Arora and Montenegro, 2011). Afforestation not only impacts the climate, but also enhances forest ecosystem services such as maintenance and enhancement of habitat provisioning and species richness (Brockerhoff et al., 2017). In recent decades, China has implemented large-scale afforestation programmes (Zhang et al., 2000), such as the Three-North Shelter Forest Program (Hu et al., 2021), the Grain for Green Programme (Xiao, 2014), and the Natural Forest Conservation Program (Huang et al., 2019). These ecological engineering programmes have been beneficial for water conservation (Liu et al., 2023), mitigating climate warming (Yu et al., 2020), increasing terrestrial carbon sequestration (Shi and Han, 2014), reducing water erosion risk (Wang et al., 2021), and alleviating dust storms (Tan and Li, 2015). These initiatives have significantly increased China's total forest cover from 8.6 % in 1949 to 24.02 % in 2022 (Zhang and Song, 2006; Fu et al., 2023; Moore et al., 2016), and they have contributed to 42 % of land greening in China during 2000–2017 (Chen et al., 2019).

In September 2020, the Chinese government declared a specific objective of achieving carbon neutrality before 2060 (Liu et al., 2022; Zhao et al., 2022). In pursuit of this goal, China is committed to expanding its forest area in the future, and new national afforestation plans have been introduced. For instance, the Action Plan for Carbon Dioxide Peaking Before 2030 (State Council of China, 2021) outlines China's target to increase forest cover to 25 % by 2030. The National Forest Management Planning (2016–2050), issued by the State Forestry Administration of China in 2016, set the afforestation target of about 73.78×10⁴ km² for the period 2020–2050 in China. Such extensive afforestation in the future would lead to land cover conversion from non-forestland to forestland, potentially causing a series of effects through the aforementioned biogeophysical and biogeochemical processes. It is important to note that the effects of afforestation are highly dependent on the afforestation location. For example, tropical afforestation may yield greater cooling effects than boreal afforestation (Arora and Montenegro, 2011). However, recent studies find that the benefits of afforestation may be overestimated, which is controversial, because the responses of the global carbon cycle to anthropogenic land use change are uncertain (Bastin et al., 2019; Veldman et al., 2019; Lewis et al., 2019). There is limited climate change net mitigation potential if tree planting occur in water-limited locations, such as in drylands (Rohatyn et al., 2022). It is therefore imperative to strategically allocate the national planned afforestation area to specific regions and project the possible land cover changes resulting from afforestation.

Various studies have examined the climatic effects of future afforestation scenarios (Abiodun et al., 2013; Naik and Abiodun, 2016; Diasso and Abiodun, 2018; Odoulami et al., 2019; Zhang et al., 2022). For example, Odoulami et al. (2019) fully replaced savanna areas (between 8 and 12° N) with evergreen broadleaf trees over West Africa to study the climate effects of future afforestation. An obvious increase in total annual precipitation was found over the afforested area. Similarly, Abiodun et al. (2013) employed random afforestation scenarios to replace 25 %–100 % of the current land cover in Nigeria and found a local cooling effect. In summary, these studies mostly employed idealistic and hypothetical afforestation scenarios and neglected the future climatological suitability of forests. In addition, process-based dynamic global vegetation models (DGVMs), such as LPJ-GUESS, have been extensively used to explore the responses of potential natural vegetation distribution to climate change (Hickler et al., 2012; Verbruggen et al., 2021), and they are also useful tools for quantifying future afforestation scenarios (Krinner et al., 2005; Horvath et al., 2021). The DGVMs driven by meteorological data generally consider complex biogeophysical, biogeochemical, and physiological processes, such as evapotranspiration, carbon–nitrogen interactions, and photosynthesis (Cramer et al., 2001). Given that the mathematical representations of these processes and their parameters as well as future meteorological scenario data from global climate models (GCMs) have large uncertainties, their overlap may yield greater uncertainties (Jiang et al., 2012; Martens et al., 2021).

The impact of future climate change is one of the biggest challenges. Previous studies (de Lima et al., 2022; Hinze et al., 2023) explored the responses of potential vegetation distribution to future climate change based on climate–vegetation models forced by the climate projection data of the GCM. However, the resolution of the raw GCM is too coarse (∼ 100–300 km) to describe the finer land surface features at a regional scale (Varney et al., 2022; Turner et al., 2023; Song and Yan, 2022; Parsons, 2020). To overcome such shortcomings, downscaling techniques are widely used to translate GCM output to regional high-resolution data. Statistical downscaling involves the establishment of statistical relationships between local climate variables and coarsely resolved atmospheric fields (Wilby and Dawson, 2013). However, it is not clear whether this historical statistical relationship is always stable in future climate scenarios. Moreover, statistical downscaling cannot ensure the physical consistency among meteorological variables. By contrast, physically based dynamical downscaling using a regional climate model (RCM) nested within a GCM could provide high-resolution climate simulations (Giorgi and Mearns, 1999; Mishra et al., 2014). Physical consistency is crucial for identifying potential afforestation regions due to the multiple meteorological variables involved. Previous studies (Liu et al., 2020a; Bowden et al., 2021) have employed the dynamical downscaling approach to quantify the climatological suitability of natural vegetation. However, previous studies (Niu et al., 2019; Wu and Gao, 2020) used the raw GCM outputs as the lateral boundary conditions (LBCs) of RCM. It is well known that raw GCM outputs have some uncertainties, and the accuracy of LBCs is the most critical factor affecting the performance of dynamical downscaling due to the underlying biases propagated into RCM through the LBCs (Sato et al., 2007; Moalafhi et al., 2017; Karypidou et al., 2023). Therefore, high-accuracy LBCs are the key to obtaining robust future potential vegetation types. Correcting the GCM outputs before dynamical downscaling is necessary for reducing the underlying uncertainty.

By taking into account the aforementioned background factors, this study aims to map the future afforestation distribution in China. It is highlighted that the results are constrained by both the national afforestation plan and future climate change. The national afforestation plan determines the total afforestation area of each province, and climate change determines where it is suitable for forest growth. Section 2 introduces the methodology. The main results are described in Sect. 3. The discussion and conclusions are presented in Sects. 4 and 5, respectively.

2 Method

2.1 Data sources

This study used three categories of data: (1) ground meteorology measurements and satellite-observed land use/land cover data, (2) national planned afforestation area data, and (3) climate modelling data from GCM as well as ERA5 reanalysis data.

2.1.1 Ground meteorology measurement data and land use/land cover data

This study used observation data of air temperature and precipitation at 2 m from the CN05.1 dataset (Wu and Gao, 2013). This dataset has a spatial resolution of 0.25°×0.25° and a temporal resolution of days from 1995 to 2014. The dataset was produced by interpolating more than 2400 meteorological stations in China using the “anomaly approach”. The CN05.1 dataset has been widely used to evaluate the performance of RCM simulations in China (Yu et al., 2015; Huang and Gao, 2018; Yan et al., 2019; Gao et al., 2023).

Land use type is a key parameter of RCMs (Mallard and Spero, 2019; Yan et al., 2021). This study used the Moderate Resolution Imaging Spectroradiometer (MODIS) land cover type dataset (MCD12Q1) for the year 2020 (Fig. S1 in the Supplement), with a spatial resolution of 500 m (Friedl et al., 2010). The MCD12Q1 features a 17-class International Geosphere-Biosphere Programme (IGBP) classification scheme (Loveland et al., 2000). It could match the first 17 default categories of land use with the Weather Research and Forecast (WRF) model (Table S1 in the Supplement). The MCD12Q1 is highly accurate globally, with an overall accuracy of approximately 75 % (Friedl et al., 2010; Sulla-Menashe et al., 2019). It was widely used to investigate land use and land cover change (You et al., 2020; Hou et al., 2022) and served as the lower boundary conditions for climate modelling (Yu et al., 2017; Ge et al., 2020; Zhao et al., 2021).

2.1.2 National planned afforestation area data

This study also used the national planned afforestation area data, which were from the National Forest Management Planning (2016–2050) (NFMP) released by the State Forestry Administration of China (2016). The NFMP presented the total national afforestation area of 73.78×10⁴ km² (equivalent to an increase in China's forest cover by 7.7 %) and the area corresponding to each province between 2020 and 2050 (Fig. 6d). The NFMP was utilized as a policy constraint to identify the future afforestation domain in China.

2.1.3 Climate modelling data and ERA5 reanalysis data

To select the optimal LBCs from GCM, Song et al. (2023) comprehensively evaluated the performance of GCM in the Coupled Model Intercomparison Project 6 (CMIP6). It was reported that the MPI-ESM1-2-HR model from the Max Planck Institute outperforms all other GCMs in East Asia. When compared with other CMIP6 models, the MPI-ESM1-2-HR model was reported to also have a higher performance in simulating various climatic variables such as the sea surface temperature (Bhattacharya et al., 2022), mean temperature (Karim et al., 2020), total precipitation (Kamruzzaman et al., 2022), and large-scale circulation (Han et al., 2022). The main configuration of the MPI-ESM1-2-HR model utilized in this study comprised the coupling of the atmospheric (ECHAM6.3) and ocean model (MPIOM version 1.6.2) and the JSBACH land surface scheme and HAMOCC ocean biogeochemistry model with a spatial resolution of 0.9375°×0.9375° on a latitude–longitude grid; the model is described in detail by Müller et al. (2018). The MPI-ESM-MR model involved in CMIP5, which was the precursor of the current MPI-ESM1-2-HR model, was widely used as the LBCs to force RCMs in carrying out finer-resolution climate simulation (Kebe et al., 2017; Ozturk et al., 2018; Crespo et al., 2023). It is well known that there are several Shared Socioeconomic Pathways (SSPs) for future climate projections in the CMIP6. Here, we used the climate projections of the MPI-ESM1-2-HR model under the middle-of-the-road pathway (i.e. the SSP2-4.5 scenario), which represents the development path most likely to occur (O'Neill et al., 2016).

The ERA5 reanalysis is the fifth-generation global reanalysis product developed by the European Centre for Medium-range Weather Forecast (ECMWF) (Hersbach et al., 2020). The state-of-the-art reanalysis data assimilated multi-source data including from ground-based meteorological measurements, satellite observations, and atmospheric sounding based on a 4D-var ensemble data assimilation system (Hersbach et al., 2020). The 6-hourly ERA5 reanalysis data with a spatial resolution of 1.0°×1.0° from 1994 to 2014 were also used as the LBCs. Climate variables for ERA5 reanalysis data and the MPI-ESM1-2-HR model include atmospheric fields (air temperature, specific humidity, zonal wind, meridional wind, geopotential height) and surface fields (i.e. sea surface temperature, surface pressure, soil temperature and moisture).

Despite performing better than other GCMs, the MPI-ESM1-2-HR model still exhibits biases. Hence, the corrections of climate mean and variance were carried out using the method reported by Xu and Yang (2012) according to Eqs. (1) and (2). The ERA5 data were used as a reference to correct the MPI-ESM1-2-HR model outputs. The MPI-ESM1-2-HR model outputs were interpolated into grid cells of 1.0°×1.0° using the bilinear interpolation method to match the ERA5 grid cells. The bias-corrected 6-hourly data of the MPI-ESM1-2-HR model retained the same means and variances as the ERA5 data (Figs. S2 and S3). This bias-corrected approach was applied to the atmospheric and surface fields:

$$\begin{array}{ll}\text{(1)}& {\displaystyle }{H}_{\text{cor}}=& {\displaystyle }{D}_{\text{GCM\_H}}\times {\displaystyle \frac{{\text{SD}}_{\text{ERA}}}{{\text{SD}}_{\text{GCM}}}}+M_{\text{ERA}}{\displaystyle }{F}_{\text{cor}}=& {\displaystyle }{D}_{\text{GCM\_F}}\times {\displaystyle \frac{{\text{SD}}_{\text{ERA}}}{{\text{SD}}_{\text{GCM}}}}+M_{\text{ERA}}\\ \text{(2)}& {\displaystyle }& {\displaystyle }+(M_{\text{GCM\_F}}-M_{\text{GCM\_H}}),\end{array}$$

where, H_cor and F_cor are bias-corrected data of 6-hourly MPI-ESM1-2-HR models over the historical period (1994–2014) and future period (2040–2060), respectively. D_{GCM_H} and D_{GCM_F} indicate anomaly by referring to the historical and future mean of MPI-ESM1-2-HR modelling, respectively. SD_ERA and SD_GCM indicate the standard deviation of ERA5 and MPI-ESM1-2-HR simulations during the historical period, respectively. ${\text{SD}}_{\text{ERA}}/{\text{SD}}_{\text{GCM}}$ denotes the variance-adjusted term. M_ERA denotes the climatological mean of ERA5 data during the historical period, and M_{GCM_F}−M_{GCM_H} indicates the mean future climate change projected by MPI-ESM1-2-HR.

2.2 Methodology

The study consists of three steps. As shown in Fig. 1, the first step is to carry out dynamical downscaling and prepare finer-resolution climate data; the second step is to run the Holdridge life zone model to identify forest-suitable lands under future climate change scenarios; finally, the third step is to allocate the national afforestation plan area to 25 km, grid cells by taking into account the climatic suitability for forests.

Figure 1 Outline for mapping the spatial distribution of future afforestation in China (PVD: potential vegetation domain; WRF: Weather Research and Forecast model; MPI: MPI-ESM1-2-HR model).

2.2.1 Dynamical downscaling of GCM outputs

In this study, the WRF model served as an RCM and was utilized to obtain high-resolution simulations (Skamarock et al., 2019). As an open-source community mesoscale numerical model, the WRF model has generally been used to investigate regional climate modelling (Wang and Kotamarthi, 2015; Cardoso et al., 2019; Moustakis et al., 2021), weather diagnosis (Ullah and Shouting, 2013; Lu et al., 2021), numerical weather prediction (Case et al., 2008; Zheng et al., 2016), and land–atmosphere interactions (Wang et al., 2013; Zhang et al., 2020, 2021). Specifically, the WRF model has been demonstrated to reproduce well the historical spatiotemporal characteristics of temperature (Politi et al., 2021), precipitation (Moustakis et al., 2022), and biome classification (Zevallos and Lavado-Casimiro, 2022), and it can successfully project the changes in temperature and precipitation over China (Hui et al., 2018). The WRF model configurations and physical parameterization (Hu et al., 2015) used in this study are detailed in Table 1. The simulation domain is shown in Fig. 2.

Table 1 Model configurations and physical parameterization for WRF simulations.

Figure 2 Model domain with topography.

The last two decades (from 1994 to 2014) were considered the historical period in this study because the historical simulation for GCM is up to 2014. Given that the NFMP is implemented for afforestation up to 2050, the simulation for the future period covers the decade around 2050, from 2040 to 2060. Three 21-year numerical experiments were performed using the WRF model (Table 2). The first two experiments, HIS_ERA and HIS_MPI, simulated the historical climate change (1994–2014) using the ERA5 analysis and MPI-ESM1-2-HR models as LBCs and default land use, respectively. The future climate change experiment (FUT_MPI) used the 2020 MCD12Q1 land cover in simulating the future period (2040–2060). All the WRF experiments were run for 21 years (1994–2014 and 2040–2060), but the first year of each experiment (1994 and 2040) as spin-up time was discarded. The remaining 20-year period (1995–2014 and 2041–2060) was analysed. We compared the HIS_MPI and HIS_ERA experiments to validate the simulation performance. The FUT_MPI experiment generated a high-resolution future climate dataset under the SSP2-4.5 scenarios.

Table 2 Detailed WRF numerical experiment design.

2.2.2 Identification of forest-suitable lands under future climate change scenario

The distribution of terrestrial ecosystems is directly affected by some main climate factors (i.e. temperature) (Piao et al., 2011; Tatli and Dalfes, 2016). Therefore, the impact of future climate change on forest-suitable lands is a further factor to be analysed. It is noted that the forest-suitable lands in this study indicate the area of the potential forestation domain (PFD). Climate–vegetation models can describe the relationship between the potential vegetation domain (PVD) and the climatic conditions (Dan et al., 2005; Kummu et al., 2021; Anwar and Diallo, 2022). Among a series of climate–vegetation models – such as the Holdridge life zone (HLZ) model (Holdridge, 1947), BIOME4 model (Kaplan, 2001), BOX model (Box, 1981), LPJ-DVGM model (Sitch et al., 2003), MAPSS model (Neilson et al., 1992), and IBIS model (Foley et al., 1996) – the HLZ model is a classification model for simulating the correlation between the potential terrestrial ecosystem types and climate change based on the combination of key climate variables (Holdridge, 1947). In recent years, the HLZ model has been globally well accepted and used to quantitatively identify the impacts of climate change on the distribution of PVD at the global (Elsen et al., 2022; Navarro et al., 2022), continental (Fan et al., 2019), and regional scales, e.g. China (Fan and Bai, 2021; Li et al., 2022). Therefore, the HLZ model was used to obtain the spatial pattern of forest-suitable lands during 2041–2060 under the SSP2-4.5 scenario over China.

The HLZ classification system requires daily temperature and monthly precipitation to obtain three bioclimatic variables: annual average biotemperature (AT), annual total precipitation (TP), and potential evapotranspiration ratio (PE). The output of the FUT_2020 experiment provides these meteorological variables. The HLZ model is estimated with a specific calculation formula as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }\text{AT}\left(t\right)={\displaystyle \frac{{\sum }_{j=\mathrm{1}}^{n}T(j,t)}{n}},\text{(4)}& {\displaystyle }\text{TP}\left(t\right)=\sum _{j=\mathrm{1}}^{n}P(j,t),\text{(5)}& {\displaystyle }\text{PE}\left(t\right)={\displaystyle \frac{\mathrm{58.93}\text{AT}\left(t\right)}{\text{TP}\left(t\right)}},\text{(6)}& {\displaystyle }\text{HLZ}\left(t\right)=\sqrt{\begin{array}{c}\left(\text{TEM}\right(t)-T_{i\mathrm{0}}{)}^{\mathrm{2}}+(\text{PER}\left(t\right)-P_{i\mathrm{0}}{)}^{\mathrm{2}}\\ +\left(\text{PET}\right(t)-E_{i\mathrm{0}}{)}^{\mathrm{2}}\end{array}},\end{array}$$

where, AT(t), TP(t), and PE(t) are the AT (°C), TP (mm), and PE for each grid in the period t, respectively. T(j,t) and P(j,t) are the mean temperature with values above 0 and below 30 °C and the total precipitation on the jth day in the period t, respectively; n is the number of days in a year. TEM(t)=ln AT(t), PER(t)=ln TP(t), PET(t)=ln PE(t); T_i0, P_i0, and E_i0 are the reference values of the classification scheme of the AT logarithm, TP logarithm, and PE logarithm, respectively, at the central point of the ith potential vegetation type in the HLZ model classification scheme. HLZ(t) is the ith potential vegetation type in the period t. A low HLZ value indicates greater potential vegetation. Fan et al. (2019) improved the HLZ model and revised the classification scheme, which was applied successfully to Eurasia. In this study, the reference values of the classification scheme were used to quantify the distribution of potential vegetation types in China (Table S3); for more details, the reader is referred to Fan et al. (2019). Compared to the actual vegetation types, the HLZ model can reproduce the potential forest distribution and grassland–forest geographical boundary well (Fig. S4).

2.2.3 Approach for allocation of new afforestation areas

In this section, we designed an approach to allocate the new afforestation area for each province into grid cells. To obtain plausible afforestation scenarios, the overall principles were that future afforestation areas should consider both future climate change and the national afforestation plan. The specific details are as follows:

• 1.

  The final total afforestation area should be consistent with the NFMP.

• 2.

  Present forestland, cropland, urban, wetland, and water body areas must not be encroached on. If the demand of the NFMP cannot be met, we consider minimizing encroachment on cropland. This can establish the concept of sustainable development as well as avoid repeated afforestation in the future (Zomer et al., 2008). The present land cover dataset for the year 2020 is based on MCD12Q1.

• 3.

  After afforestation, China's cultivated land area is not expected to fall below 121.67×10⁴ km² according to the requirements of the National Land Planning Outline (2016–2030) (State Council of China, 2017). This ensures that the cultivated land area stays within the “red line” and enhances people's welfare.

• 4.

  Afforestation is implemented in areas where the potential vegetation type is forestland in the context of future climate, according to the output of the HLZ model. This measure could ensure that future climate conditions are suitable for the growth of forests.

• 5.

  Areas with a low HLZ value are given priority afforestation. The HLZ metric is a comprehensive metric considering the biotemperature, precipitation, and potential evapotranspiration ratio. A low HLZ value means a greater chance of the area being potential forestland according to the HLZ model.

3 Results

3.1 Model evaluation

We evaluate the performance of the key climate variables and PVD based on the observations and two WRF simulations (HIS_ERA and HIS_MPI). The performance of the WRF simulation is quantified by the bias, mean absolute error (MAE), and spatial correlation coefficient (R) for the bioclimatic variables (AT, TP, and PE). Larger R values and smaller bias and MAE values indicate better performance. Figure 3 illustrates the spatial patterns of the observed and simulated multi-annual averaged (1995–2014) AT, TP, PE, and HLZ types in China. The WRF simulation can reproduce well the spatial distribution of the observations with an increasing northwest to southeast temperature and precipitation gradient. However, the underlying bias still remains against the observations (Figs. S5–S7). A more detailed inspection of the scatterplots shows that the spatial correlation coefficient between the observation and simulations (HIS_ERA) is 0.982 for AT, 0.795 for TP, and 0.754 for PE (Fig. 4). The simulated AT is generally underestimated in most regions, with the national-average bias of −0.974 °C. Consistent with the previous studies (Meng et al., 2018), the largest cold biases are in the Tibetan Plateau and complex terrain regions (Fig. S5), with a bias of more than −3.6 °C, which could be attributed to the poor simulation of the snow–ice albedo feedback process (Ji and Kang, 2013). The simulated AT is relatively better in eastern China. The WRF simulation generally overestimates TP in most regions with a national-average bias of 92.883 mm (Fig. 4d). The wet bias could be attributed to inappropriate parameterization schemes (Ou et al., 2020; Zhao et al., 2023), coarse horizontal resolution (Lin et al., 2018; Rahimi et al., 2019), and inappropriate land surface processes associated with soil moisture and freezing–thawing (Fu et al., 2020; Yang et al., 2018). However, the scatterplot dispersion displays that the simulated TP exceeding 1200 mm in southern China is underestimated (Fig. S5). It is not surprising that the temperature is well modelled, but the simulation capacity of precipitation-related variables is modest for the WRF model (Gao, 2020). It should be noted that the HIS_ERA simulation exhibits a highly consistent representation to that of HIS_MPI. The cross-correlations for three climate variables between the HIS_ERA and HIS_MPI simulation show a high spatial correlation coefficient, and the scatter distribution is very close to the 1:1 line (Fig. 4).

Figure 3 Spatial pattern of annual average biotemperature (AT), annual total precipitation (TP), potential evapotranspiration ratio (PE), and potential vegetation domain from the Holdridge life zone (HLZ) model based on observations (a, d, g, j), HIS_ERA simulations (b, e, h, k), and HIS_MPI simulations (c, f, i, l) during the period 1995–2014.

Figure 4 Scatterplots of the annual average biotemperature (AT), annual total precipitation (TP), and potential evapotranspiration ratio (PE) for each grid against the observations and HIS_ERA, the observations and HIS_MPI, and HIS_MPI and HIS_ERA. HIS_ERA and HIS_MPI indicate the WRF simulation driven by the ERA5 reanalysis data and the bias-corrected MPI-ESM1-2-HR model, respectively. The observation is derived from the CN05.1 dataset. Evaluation indexes include the bias, mean absolute error (MAE), and spatial correlation coefficient (R). The dotted black line indicates the 1:1 line.

The observed and simulated results of PVD are shown in Fig. 3j–l. The kappa statistic is applied to validate the observed and simulated accuracy of the PVD map from the HLZ model (Cohen, 1960). The kappa coefficient ranges from 0 to 1.0, and the degree of agreement differs across these ranges. According to the description of Landis and Koch (1977), a kappa coefficient range of 0–0.2 is considered slight agreement, 0.21–0.40 is fair agreement, 0.41–0.60 is moderate agreement, 0.61–0.80 is substantial agreement, and 0.81–1.00 is almost perfect agreement. Overall, the WRF_ERA simulation could reproduce the distribution of PVD well in China. However, some minor differences in vegetation types are found. For example, in the northeast region of China, the WRF simulation could not precisely reproduce the observed extent of steppe types. Such misclassified zones could be due to the fact that the model overestimated the precipitation exceeding 220 mm in the transition zone of a dry–wet climate (Fig. S5d); thus, the vegetation types are changed from steppe to cool temperature forest. Other disagreement types are found in southern China, where the observed subtropical forest expands northward up to 32° N. However, the simulation results reduce the extent. The dry bias of precipitation simulation in southern China could explain the source of uncertainty. Despite the limited ability displayed by the models, the overall accuracy based on the kappa coefficient indicates a substantial agreement between the observation and the WRF simulation. When compared with the observations, the kappa coefficient is 0.648 for HIS_ERA and 0.662 for HIS_MPI. This suggests a perfect agreement (kappa coefficient = 0.962) of the PVD between the HIS_MPI and HIS_ERA (as shown in Fig. 3k and l). The implementation suggests that the bias-corrected MPI-ESM1-2-HR model can replace the ERA5 reanalysis data as the LBCs of WRF to obtain similar accuracy in high-resolution simulations. Therefore, in the following analysis, the WRF model forced by the bias-corrected MPI-ESM1-2-HR model will be applied to project future climate change.

3.2 Future potential vegetation cover

For the future simulation, the three key variables (AT, TP, and PE) of the HLZ model are obtained from the FUT_MPI experiment. The projected spatial distribution of PVD is presented in Fig. 5a. The most dominant vegetation types are forest, polar/tundra, and desert/scrub, accounting for 57.1 %, 20.1 %, and 17.7 % of the total area of China, respectively. The forest types are located in eastern China, characterized by a latitudinal distribution. The potential forest types from north to south are mainly cool temperate forests, warm temperate forests, and subtropical forests, in that order. This may be explained by the fact that temperature is the critical factor in defining the forest types due to the sufficient precipitation in eastern China.

Figure 5 Projected spatial pattern of (a) potential vegetation types from the HLZ model under the SSP2-4.5 scenario in the future period (2041–2060) from the FUT_MPI simulation and (b) area changes across historical baseline (1995–2014) and future period, where the calculations are based on the FUT_MPI simulation versus the HIS_MPI simulation.

Flow diagrams are useful tools for depicting precise changes in vegetation types, displaying whether the vegetation types are shifting and in which direction. The projected changes in the vegetation types for the area are shown in Fig. 5b. The results indicate that under future climate change, the PVD changes correspondingly. A total area of 10.4 % will be shifted in China. The northward expansion of subtropical forests replaces warm temperate forests, with an area of approximately 30.6×10⁴ km², which is considered the largest shift type (Fig. S8). In addition, projected future increases in temperature and precipitation have caused some non-forestland areas to transition into forested lands (Fig. S9). For example, in western China, areas that are polar/tundra and steppe at present have transitioned into cold temperate forest and cool temperate forest in the future, respectively, and the shifted area is 18.4×10⁴ and 1.7×10⁴ km², respectively.

Overall, the PFD (1995–2014) covers approximately 500.75×10⁴ km². It is projected to expand to 518.25×10⁴ km², experiencing an increase of around 17.5×10⁴ km² (about 3.49 %) within the period 2041–2060 under the SSP2-4.5 scenario. In eastern China, the main transition is interconversion between forest types. In western China, some non-forestland types turn into forest types. These changes indicate that the forest-suitable region would be modified under future climate change, and it is necessary to consider the climatic contexts in terms of future large-scale afforestation.

3.3 Identification of future potential afforestation location in China

According to the approach of the new afforestation allocation in Sect. 2.2.3, we mapped the future afforestation distribution of China. First, historical open-space regions for afforestation were identified. We excluded some ineligible regions, including present forestland, cropland, urban, wetland, and water body areas based on the 2020 MCD12Q1 land cover data (Fig. S1); the remaining regions were considered open-space regions for afforestation (Rohatyn et al., 2022). The results show that the total area of open-space regions is about 612.88×10⁴ km² in China, with the majority located in southern and western China (Fig. 6a).

Figure 6 Spatial distribution of (a) historical open-space region for afforestation, (b) future potential forestation domain (PFD) from the HLZ model considered to be forest-suitable lands, (c) potential afforestation region constrained by climate change, (d) national planned afforestation area in the individual provinces from the NFMP, and (e) Chinese vegetation regionalization map.

The second step was to determine the distribution of future PFD. We used the map of potential vegetation derived from the outputs of the HLZ model (Fig. 5a) to select the forest type grids as future PFD under the SSP2-4.5 scenarios during 2041–2060. The future PFD was considered to be forest-suitable lands constrained by future climate conditions. The forest-suitable lands are mainly located in eastern China (Fig. 6b). The corresponding annual total precipitation is over 353.6 mm among the selected grids.

Then, we combined the historical open-space region (Fig. 6a) with the future PFD (Fig. 6b), which enabled us to obtain the future potential afforestation areas (Fig. 6c). These regions provide suitable climate conditions for forest growth and can be utilized for afforestation implementation in the context of future climate change. The total area of potential afforestation is approximately 191.33×10⁴ km².

There is no doubt that the potential afforestation area is extensive and unrealistic. Thus, according to the national tree planning policy, we further restricted the afforestation area. The NFMP released by the State Forestry Administration of China (2016) included the total area of planning afforestation in each province during 2020–2050 (Fig. 6d) and was considered a reference for future afforestation design. It notes that the potential afforestation area for individual provinces is usually larger than the national planned afforestation area (Table S2). Thus, we further constrained the potential afforestation areas following the HLZ value. Specifically, we sorted the HLZ value for each province in the potential afforestation region in ascending order (Fig. 6c); low HLZ values were given priority afforestation. We calculated the total afforestation area sequentially grid by grid, until it satisfied the NFMP policy requirements. The approach of the total afforestation area for each province is calculated based on Eq. (7):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\displaystyle }\text{Area}=(& {\displaystyle }\mathrm{0.55}{N}_{\text{Woody savannas}}+\mathrm{0.80}{N}_{\text{Savannas}}\\ {\displaystyle }& {\displaystyle }+N_{\text{Grasslands and croplands}})\times r^{\mathrm{2}},\end{array}\end{array}$$

where Area indicates the total afforestation area; r indicates the spatial resolution (here, r equals 25 km); and N indicates the number of afforestation grids in the historical land cover. The land cover types represent the area used for afforestation. Given the tree cover for woody savannas and savannas is 30 %–60 % and 10 %–30 % according to the IGBP classification scheme (Table S1), it means that approximately 45 % and 20 % of the grid area for woody savannas and savannas, respectively, has already been covered by forests in the historical period. Thus, to avoid repeated afforestation, the coefficients 0.55 and 0.80 were set.

It is worth noting that the planned afforestation area is larger than the potential afforestation area in Henan and Shandong provinces (Table S2). A small amount of cropland had been scheduled for afforestation to meet the national afforestation demand. The occupied croplands are mainly located in mountain areas, where the regions are highly suitable for forest growth. With such an afforestation scenario design, 125.33×10⁴ km² croplands in China are still available for cultivation. It is also away from the protection “red line” of 121.67×10⁴ km², released by the National Land Planning Outline (2016–2030) (State Council of China, 2017).

A Chinese vegetation regionalization map (Wu et al., 1980) was used to identify the forest types within each grid (Fig. 6e). Finally, the distribution of future potential afforestation regions in China is shown in Fig. 7. The findings show that the probable locations for future potential afforestation areas in China are around and to the east of the Hu Line (a geographical division of climate zone, population density, and economic development in China, stretching from Heihe to Tengchong). Due to afforestation, the land cover would be modified. In northern China, the main conversion types are from grasslands to deciduous broadleaf forests, and the region has the largest land conversions in China, accounting for 40 % of the new afforestation area. The most intensive conversion provinces are Shanxi and Shaanxi. In southwest China, the dominant conversions are from woody savannas and savannas to evergreen broadleaf forests. These conversions account for 26 % and 16 % of the new afforestation area, respectively. These land use conversions are mainly located in southwest China, such as Yunnan province, Sichuan province, and Guizhou province. Overall, the final total afforestation area in China is approximately 73.64×10⁴ km², consistent with the NFMP (73.78×10⁴ km²). Therefore, for each province within the future afforestation region, we applied the aforementioned approach to ensure that the total afforestation area of individual provinces and the extent were consistent with the national policies and future climate conditions, respectively.

Figure 7 Map of future potential afforestation distribution and land shift types constrained by both the national afforestation plan and climate change. Forest types from IGBP include evergreen needleleaf forests (ENF), evergreen broadleaf forests (EBF), deciduous needleleaf forests (DNF), deciduous broadleaf forests (DBF), and mixed forests (MF). The dotted black line indicates the Hu Line.

4 Discussion

The most probable geographical distribution of future potential afforestation regions in China was investigated in this study. The total afforestation area in this study (73.64×10⁴ km²) is larger than that in existing studies. For example, Zhang et al. (2022) reported an obvious increase in potential forestation lands by 33.1×10⁴ km² under future climate scenarios (2070s) using a machine learning approach to predict the ecological niche of forests. Xu (2023) reported that the area of prioritized potential forestation land was about 66.61×10⁴ km² in 2020 through spatial overlay analysis of multiple factors (i.e. climate, transportation, topography, land use). However, the effects of future climate change and the national afforestation plan are ignored. Our results show that forest-suitable lands will increase by 17.5×10⁴ km² under the SSP2-4.5 scenario compared with the historical period. The dataset would be valuable for studying the effects of future afforestation on carbon budget, ecosystem services, water resources, and surface climate.

Our findings indicate that future afforestation in China would mostly be located around and to the east of the Hu Line, consistent with Zhang et al. (2022). The area near the Hu Line is a transition zone characterized by dry–wet, agro-pastoral, and grassland–forest areas. This transition zone is highly sensitive to climate change (Li et al., 2015). Due to moisture limitations, the historical forest distribution is mainly located east of the Hu Line. Crossing the Hu Line is challenging for forests (Liu, 2019). However, under future climate change, the projected results show that the temperature and precipitation in China will increase by the middle of the 21st century under the SSP2-4.5 scenario compared with the historical period (Yang et al., 2021). A similar conclusion is also derived from our projection (Fig. S9). The response of PFD to future climate change could be slightly modified. Therefore, only a small proportion of future potential afforestation areas are in the western region of the Hu Line, such as the Loess Plateau region. This shows that afforestation planning should consider vegetation responses to future climate change.

Afforestation can provide temperature benefits (e.g. cooling the land surface), according to previous studies (Peng et al., 2014; Yu et al., 2020; Breil et al., 2024); however, the biophysical response of afforestation to temperature varies spatially. At a global scale, it is logical that afforestation causes a warming effect in high-latitude regions due to the albedo-dominant radiation effect, while it causes a cooling effect in low-latitude regions due to the evapotranspiration-dominant non-radiation effect (Bonan, 2008; Arora and Montenegro, 2011). Thus, afforestation-induced regional temperature changes depend on the net effects. Afforestation can also cause daytime cooling but nighttime warming (Yuan et al., 2022), and it can increase the surface temperature in winter but decrease it in the other seasons (Ma et al., 2017). Differential seasonal and daily responses lead to larger uncertainties in the net effects induced by afforestation. Therefore, a more realistic afforestation scenario is necessary in order to quantify the effects of afforestation on temperature under future climate change and to develop climate change mitigation policies.

Although the resolution of our dynamical downscaled simulation (25 km) is finer than raw GCM (∼ 100 km), it is difficult to meet the needs of afforestation planning in areas with complex topography. Convection-permitting climate modelling at the kilometre scale has recently been developed to reproduce better mesoscale atmospheric processes (Prein et al., 2015; Lucas-Picher et al., 2021) and obviously improve the WRF simulation, especially precipitation (Knist et al., 2020). However, improving the resolution involves higher computational costs. By contrast, statistical downscaling methods are also known to obtain high-resolution climate data with fewer computational resources (Tang et al., 2016). These methods assume that the historical relationship between local climate variables and the large-scale circulation remains fixed in the future term (Wilby and Dawson, 2013). The multi-model ensemble from the CMIP6 statistical downscaling can significantly reduce the biases compared with individual models (Gebrechorkos et al., 2019). Thus, some statistical downscaled CMIP6 datasets (Gebrechorkos et al., 2023; Lin et al., 2023; Thrasher et al., 2022), with a resolution of 0.1–0.25° covering the global area, can be applied to explore the future global potential afforestation area in later work. However, it is noted that the statistical downscaling data may have a limitation, as the covariance among the variables may not align with physical laws.

This study may have some limitations and uncertainties. Following the approach of existing studies (Ma et al., 2023; Qiu et al., 2022), we also utilized the bias-corrected LBCs in dynamical downscaling. However, the model uncertainty in the future climate projection is difficult to quantify because one GCM is nested within the WRF model. The projected result generally exhibits variations based on the choice of driving GCMs (Gao et al., 2022). This divergence can be attributed to the inherent configurations and physical parameterization of the GCM, distinct radiative forcing scenarios, and varying equilibrium climate sensitivities found in CMIP6 models (Zuo et al., 2023; Bukovsky and Mearns, 2020). For instance, the high-emission scenario could lead to higher temperature and stronger precipitation in China relative to the middle-emission scenario (Yang et al., 2021). Obvious differences are found in northern China, suggesting that there are greater opportunities for afforestation in semi-arid areas. Thus, the suitability of future forest lands depends on emission scenarios (Liu et al., 2020b; Elsen et al., 2022). Exploring the impacts of different SSPs on the distribution of potential afforestation regions would be an intriguing avenue for future research. To address the concerns about model uncertainty, exploring WRF forced by multiple bias-corrected CMIP6 models can help uncover the source of uncertainty. Utilizing ensemble means for downscaled climate simulation would contribute to a more robust projection. Additionally, the selection of different physical parameterization schemes in the WRF model can also influence the simulation performance (Gbode et al., 2019). Selecting the most suitable combination is beneficial for reducing the underlying bias. Out of all the factors limiting afforestation allocation, we used the HLZ value to constrain the afforestation distribution. Previous studies found that precipitation was a key meteorological factor that restricts forest distribution, especially in the mid-latitude regions (Hansen et al., 2005; Fang et al., 2005). When the areas with high precipitation were given priority afforestation, we obtained a similar future potential afforestation distribution (Fig. S10). Future studies should comprehensively consider additional factors, such as local economic development, soil physicochemical properties, and provincial tree planning policy.

5 Conclusions

This study evaluated the performance of the WRF model in simulating the PVD from the HLZ model in China during the historical period (1995–2014). The projected shifts in the potential vegetation types were explored under the SSP2-4.5 scenario during the future period (2041–2060) relative to the historical period. Based on these data, the most probable distribution of future potential afforestation was obtained by constraining both future climate contexts and national afforestation plans in China. We could draw the following main conclusions:

The output of the WRF model forced by the ERA5 analysis and bias-corrected MPI-ESM1-2-HR model could capture the spatial distribution of the PVD from the HLZ model over China through comparisons with the CN05.1 dataset during the historical period. However, the WRF simulation did not precisely reproduce the observed extent of steppe types in northeast China and of subtropical forests in southern China. Such misclassifications might be attributed to the bias of the precipitation simulation. Overall, in terms of the nationwide potential forestation domain, the WRF model could reproduce the spatial distribution well over China.

Under the SSP2-4.5 scenario, the PVD would obviously shift during 2041–2060 compared to the historical period. The largest shift type was from warm temperate forests to subtropical forests over southern China. The new forest-suitable lands would increase by about 17.5×10⁴ km² in China due to projected increases in temperature and precipitation. In addition, considering both future climate change and the national tree planning policy, we found that the probable locations for future afforestation were around and to the east of the Hu Line, with a total area of approximately 73.64×10⁴ km². The main shift types were from grasslands to deciduous broadleaf forests in northern China and from woody savannas and savannas to evergreen broadleaf forests in southwest China. The findings of this study could provide a dataset for exploring the effects of future afforestation, and this method can guide designs of future gridded afforestation regions for other countries.