The study area is the entire subtropical region of China, which is located to the south of the Qinling Mountains and Huaihe River, north of Leizhou Peninsula, and east of the Hengduan Mountains (22° -34°N, 98° -122°E). The terrain is low in the west and high in the east. China’s subtropics belong to the east coast humid monsoon area, and the region is the warmest and hottest compared to the same latitude, except for desert areas. Moreover, the rainfall in the study region is far more abundant than that in the same latitude worldwide (Yang et al., 2006). The average annual temperature ranges from -1 to 24°C, and the average annual precipitation ranges from 450 to 2125mm. China has preserved the best subtropical evergreen forest ecosystem which is the main component of China’s subtropical forests (Lin et al., 2020) and accounts for approximately 25% of China’s land area. Dominant families in such ecosystems are the Cyclobalanopsis Oerst, Castanopsis Spach, and Lithocarpus of Fagaceae.

In this study, the SPEI were calculated based on soil heat flux, meteorological data, saturated water vapor pressure, sunshine duration, wind speed, and latitude. The calculation process of the SPEI based on this algorithm is as follows (Zotarelli et al., 2010).

1). The core of the SPEI algorithm is the calculation for PET. The Penman-Monteith (P-M) formula considered heat and aerodynamics factors, and is more consistent with the measured evapotranspiration (Penman, 1948; Jensen et al., 1990). The P-M formula is as follows:

where R_s is the daily average net radiation flux of the plant surface (MJ·m^-2·day^-1); G is the soil heat flux (MJ·m^-2·day^-1); γ is the humidity constant (kPa·°C^-1); T is the average temperature(°C); U₂ is the average wind speed at 2 m per day (m·day^-1); and e_s, e_a and Δ are the saturated water vapor pressure (kPa), actual water vapor pressure (kPa) and slope of the saturated water vapor pressure-temperature curve, respectively. R_s, G, γ, and other parameters were calculated using the FAO-56 P-M updating equation recommended by the FAO (Allen et al., 1998).

2). Calculate the measurement of deficit water balance (D_i):

3). Create a series of cumulative water deficits at different time scales D n k :

where k is the time scale (generally on month) and n is the number of calculations.

4). Choose the log-logistic probability distribution (negative values can be interpreted and modeled with different shapes for the frequencies of the D series at different timescales) to normalize the D series:

In formula (4) the following parameters are defined:

where Γ is a factorial function and w₀, w₁, and w₂ are the probability-weighted moments of D_i.

where N is the number of months involved in the calculation.

5). Standardize the cumulative probability density to obtain the SPEI (the standard value of SPEI is 0, and the standard deviation is 1):

First, calculate the size of P:

①

W h e n   P ≤ 0.5 , w = - 2 l n ( P ) ,

②

W h e n P > 0.5 , P = 1 − P , w = - 2 l n ( P ) ,

In the above equations, c_{0 =} 2.515517, c_{1 =} 0.802853, c_{2 =} 0.010328, d_{1 =} 1.432788, d_{2 =} 0.189269, and d_{3 =} 0.001308.

The SPEI usually has a variety of time scales, such as monthly, seasonal, and annual. Compared to the SPEI with shorter time scales for targeting meteorological and agricultural drought, the SPEI with a longer time scale is more sensitive to hydrological drought (Ivits et al., 2014), moreover, vegetation germination in biological communities is mainly affected by accumulated precipitation in the previous 12 months (Vicente-Serrano, 2006). Therefore, in this study, the SPEI value at the 12-month scale (hereafter referred to as SPEI₁₂) was used to characterize the interannual drought in subtropical regions, and the linear regression equation was used to evaluate the trend and characteristics of SPEI₁₂ over time based on the following formula:

where slope is the changing trend (when slope >0, it indicates that SPEI₁₂ is increasing; when slope<0, then the SPEI₁₂ is decreasing), y is the number of drought years, i is the number of years (i=1, 2… n), and SNP_i is the SPEI₁₂ value of the i_th year.

On the spatial scale, the drought classification criteria ( Table 1 ) (Center et al., 2017) were used to define the drought conditions among different regions in subtropical China. For the time series, the overall mean value could not meet the criteria for drought classification due to the spatial variation of SPEI₁₂. Therefore, in this study, the threshold value of the SPEI₁₂ percentile was used to define the annual drought conditions with reference to the precipitation percentile threshold method (Zhai and Pan, 2003), as shown in Table 2 .

In this study, the carbon storage of vegetation in two types of forests in subtropical China from 1985 to 2019 was simulated with a spatial resolution of 1km and temporal resolution on a daily time steps.

The Biome-BGC model simulates the physiological and ecological processes of vegetation that control the material cycle and energy flow of the ecosystem, including canopy radiation, photosynthesis, stomatal conductance, autotrophic respiration, heterotrophic respiration, phenological dynamics, and evapotranspiration. The model can simulate the energy and carbon-nitrogen water cycles between the atmosphere, vegetation, and soil of the terrestrial ecosystem in daily steps to estimate the storage and flux fluxes among carbon, nitrogen, and water pools (Running, 1993).

The input data to the Biome-BGC model included vegetation abundance data, topographic data, soil data, and meteorological data, which has been described in Section 2.2.

This study used 67 physiological and ecological parameters to run the model. At present, few studies have focused on the physiology and ecology parameters of vegetation. In this study, the proportion of nitrogen in the Rubisco enzyme, litter coefficient at the leaf replacement period, and the carbon-nitrogen distribution ratio of each part were obtained by an iterative method (Lu et al., 2016). In addition, the default ENF or EBF parameters provided by White (White et al., 2000) were adopted if the parameters could not be determined through the literature. Some parameter values are shown in Table 3 .

First, the model was spun-up, meaning that the carbon, nitrogen, and water storage of the ecosystem when the annual change in the soil carbon pool was less than 0.0005 kg C·m^-2 (Thornton, 2010) was used as the initial condition of the simulation, and then the vegetation carbon storage of the EBF and ENF in subtropical China from 1985 to 2019 was simulated.

Based on the simulation results, the linear regression analysis was used to calculate the trend of simulated carbon storage. The equation is similar to Eq. 13, where slope is the changing trend (when slope >0, carbon storage is increasing, and when slope<0, carbon storage is decreasing), y is the number of simulated years, i is the number of years (i=1, 2,…, n), and SNP_i is the value of carbon storage in the i _th year.

Forest ecosystem resistance (Rt) represents the ability of a forest ecosystem to maintain its original state under drought disturbance, whereas forest ecosystem resilience (Rs) describes the ability of a forest to recover from drought disturbance to a normal state (Van Ruijven and Berendse, 2010; De Keersmaecker et al., 2014). These parameters represent functions of the stability of an ecosystem (De Keersmaecker et al., 2014). In this study, formulas 14 and 15 were used to quantitatively explain the Rt and Rs of subtropical forests to drought, where the values of Rt and Rs are unitless to facilitate the comparison of the stability of forest ecosystems with two different levels of productivity.

where Y _n is the carbon storage in the normal year from 1985 to 2019, Y _e is the carbon storage in the extreme drought year, and Y _{e + i} is the forest carbon storage in year i after the event. Since the effects of drought on forest growth may last for several years and lead to legacy effects (Anderegg et al., 2015; Schwalm et al., 2017), i =1, 2, and 4 was used in this study to quantitatively analyze the resilience of carbon storage in the first, second, and fourth years after drought (in addition, the water balance of the forest in the four years after drought should also be considered). A higher Rt value indicates that the forest is more resistant to drought and a higher Rs value means stronger forest resilience (Huang and Xia, 2019).

In this study, we first analyzed the spatiotemporal distribution trend of the SPEI₁₂ in the subtropical region of China from 1952 to 2019. Extreme drought years were determined based on the drought time series. Combined with the simulation results of carbon storage in subtropical forests, the resistance and resilience of forests affected by extreme drought were analyzed. As the frequency of drought changes, the acclimation of forests to drought may also change accordingly (Anderegg et al., 2015; Isbell et al., 2015; Anderegg et al., 2020). Therefore, whether continuous drought occurred before and after the drought year was determined in this study, and whether continuous drought affects forest resistance and resilience is also a significant step in this research and discussion.

The time series of SPEI₁₂ in the subtropical region from 1952 to 2019 are shown in Figure 3A . During the 68 years, subtropical forests in China experienced five extreme drought events, which occurred every 13.6 years. According to the spatial variation trends in drought ( Figure 3B ), SPEI₁₂ showed a decreasing trend in 62.13% of the region, with an overall decrease of 0.036(10 a)^-1. Figure 3B showed the trends and characteristics of SPEI₁₂ with time using the linear regression equation (Eq. 13). The regions with a serious downward trend were mainly distributed in most of the western subtropical region, especially in parts of the Sichuan and Yunnan provinces, with a downward trend of more than 0.15 (10a)^-1.

According to Figure 3A , two extreme drought events occurred after 1985 in 2009 and 2011. Because the SPEI₁₂ value was lower in 2011, 2011 was defined as the extreme drought year in this study. Therefore, this study focused on the resistance of subtropical forests in 2011 to drought and their resilience at the first to fourth year after the 2011 extreme drought, which covers the period from 2011 to 2015.

To compare the different drought conditions of the two forest types, we calculated the pixel proportions of the two forest types under different drought conditions ( Figure 4 ). According in Figure 4A , the areas of extreme and severe drought were concentrated in the central to southwestern region of subtropical China, where the forest cover rate was high. The results for EBF showed that 19.3% suffered from extreme drought (SPEI₁₂≤-2) and 26.6% suffered from severe drought (-2< SPEI₁₂≤-1.5), while those for ENF showed that 31.9% suffered from extreme drought (SPEI₁₂≤-2) and 21.2% suffered from severe drought condition (-2< SPEI₁₂≤-1.5). The results showed that most EBF and ENF in subtropical regions of China suffered from severe to extreme drought in 2011.

The value of SPEI12 in 2012 is -0.028, and the water balance condition is normal as shown in Figure 3 . As shown in Figure 4B , only 8.84% of EBF pixels and 1.7% of ENF pixels are under severe to extreme drought conditions in 2012, and the majority of them are under normal water balance condition.

The SPEI₁₂ value in 2013 was -0.526, which represents a moderate drought according to the time series. However, according to the spatial distribution ( Figure 4C ), only a few forest pixels in the subtropical region of China were affected by extreme to severe drought in 2013, and the subtropical forests as a whole had normal moisture status.

The SPEI₁₂ value was 0.587 in 2015, indicating the absence of drought. Further analysis showed that the vast majority of forest pixels in the EBF (95.34%) and ENF (79.56%) were in the normal to wet state without drought in 2015 ( Figure 4D ).

The spatial distribution of average vegetation carbon storage of EBF and ENF during 1985-2019 were shown in Figure 5 . As shown in the figure, the average vegetation carbon storage of EBF and ENF ranged in 124.03 - 143.45 t hm^-2 (130.58 ± 10.02 t hm^-2) and 77.21 - 82.59 t hm^-2 (78.49 ± 8.49 t hm^-2), respectively. In terms of spatial distribution, the regions with higher carbon storage of EBF were mainly concentrated in central and southeastern subtropics, where those of ENF mainly distributed in the western subtropic, such as Tibet, Yunnan and Sichuan province.

According to the analysis of spatiotemporal characteristics of subtropical drought, follows we mainly focus on the vegetation carbon storage of two types of forest during 2011-2015 for better understanding the resistance and resilience of subtropical forests to drought. The monthly time series and spatial trends of carbon storage of two subtropical forest are shown in Figure 6 . As shown in Figure 6A , the carbon storage of EBF fluctuates within a certain range, whereas the ENF fluctuates relatively little and continues to increase each year. From 2011 to 2015, the mean values of the overall annual trend of carbon storage for the two forest species in subtropical China were -0.042 (EBF) and -0.013 (ENF). The overall carbon storage of subtropical evergreen forests showed a slightly decreasing trend, and the annual variation ranged mainly between -0.25 and 0.25. In subtropical China, 60.5% of the EBF and 51.18% of the ENF pixels showed a downward trend from 2011 to 2015. Figures 6B, C shows the trends and characteristics of carbon storage changes over time evaluated using the linear regression equation (Eq. 13). The spatial distribution of the changing trend in carbon storage in the two forests ( Figures 6B, C ) was relatively complex: the two forests showed an overall downward trend; the pixels in the middle and southeast of the EBF showed an overall upward trend, and the pixels from northwest to southeast of the ENF showed an upward trend.

To evaluate the stability of forests in subtropical China to drought, we first explored the resistance of the two forest types to drought in 2011 and the resilience of the forests at 1, 2, and 4 years after 2011.

As mentioned above, changes in drought frequency affect forest resistance and resilience; therefore, this section examines the effects of continuous drought on forest resistance and resilience.

In this study, the year of extreme drought was 2011. Pixels that were or were not affected by drought in 2009, 2010, and 2011 were divided into two categories, and pixels that presented the two different situations in the three years were combined. The forest pixels that experienced three consecutive drought years (2009-2011), two consecutive drought years (2010-2011), and one drought year (2011) were screened and divided into three climatic combinations. Figure 10 shows the resistance and resilience of the two forests under the three drought combinations.

Figure 10 shows that the EBF generally has higher resistance under three different continuous drought conditions, that is, subtropical EBF in China can quickly adapt to drought. In terms of resilience, under the condition of a different combination of drought, two kinds of resilience showed similar levels of forest, and in the fourth year after the drought, two kinds of resilience of forest showed very similar levels, suggesting that the two kinds of forest will recover to normal levels after four years of drought. These results support the idea that forests need at least one year to recover from interannual droughts (Schwalm et al., 2017).

The Thornthwaite (TW) (Thornthwaite, 1948) and P-M (Monteith, 1965) formulas were commonly used to calculate the PET in the SPEI estimation. The SPEI calculated based on TW formula was easily to implement, while it may indicate excessive dry conditions due to the influence of temperature, under the significantly increased temperatures in recent years (Zhang et al., 2017; Nouri and Homaee, 2020; Guo et al., 2022). Compared to the TW formula, which only considers temperature, the P-M formula considered variety of meteorological factors, which are more complex, and is generally more consistent with actual evapotranspiration (Jensen et al., 1990), specially it enable to describe the regions influenced by aerodynamic factors (Liu and Jiang, 2015). Consider, these two methods were widely used in the drought detection researches, further validation is required to determine which algorithm is better characterized in subtropical China.

The SPEI results were limited by meteorological stations, whereas some areas had few meteorological stations. Therefore, the accuracy of the SPEI in characterizing drought conditions in the study area needs to be confirmed by further comparative analysis with more observational data. The results of the temporal analysis indicate that the SPEI₁₂ shows a significant downward trend from the 1950s to the 2010s, the interannual drought in the subtropical region of China shows an increasing trend from the 1950s to the 2010s, which is consistent with previous studies (Liu and Jiang, 2015). Spatial analysis showed a clear trend of drought in the subtropical southwest, especially in the Sichuan Basin, which is consistent with previous research findings (Li et al., 2012). A comparison of historical drought events showed that in 1982, Ningbo, Jinhua, Wenzhou, Jiujiang, and Nanchang (Zhang and Liu, 1993; Zhang et al., 2003); in 1986, Nanchang, Shangrao, Ji ‘an, and Ganzhou (Zhang and Liu, 1993; Zhang et al., 2003); in 1991, south Guangdong and Guangxi (Wang and Chen, 2012); in 1994, Suzhou and Shanghai (Zhang et al., 2003); and in 2000, moderate to severe drought in Yangzhou, Nanjing, Hefei and Anqing (Zhang et al., 2003). The simulation results of SPEI in this study were consistent with typical historical drought events ( Table 4 ).

The uncertainty of the Biome-BGC model may be divided into the uncertainty of the model structure, input variables, and uncertainty of model parameters (Li and Sun, 2018). First, the uncertainty of the model structure may be due to inadequate simulation of carbon, nitrogen, and water cycles in the ecosystem, which may lead to the difference between the simulation and observations (Churkina et al., 2003; Hidy et al., 2012; Smith and Dukes, 2013). Impacts like human activities (Du et al., 2021) (e.g., management practices and forest wildfires) on the carbon and nitrogen water cycle were not considered in this simulation, which may create uncertainty. Second, the uncertainty of the input variables may be due to errors or inadequacies in the collection and statistics of input data (Jung et al., 2007; Eastaugh et al., 2011). Although the forest abundance data used in this simulation are from fine-scale classification products with a resolution of 30 m, resampling to 1 km does not avoid pixel mixing and thus errors (Su et al., 2022). Third, the uncertainty of the model parameters may be due to their different effects of the model parameters on the output results under different conditions (Tatarinov and Cienciala, 2006; Kang, 2016; Raj et al., 2018). For some difficult-to-obtain parameters, this study was obtained by reviewing the literature and directly using the model defaults, which may cause uncertainty in the results.

However, by comparing the simulated values with the observed values in the Zhejiang forest inventory sample plots, and using evaluation indicators such as correlation coefficient and root mean square error to analyze the simulation accuracy, we found that the simulated carbon storage values of both forest types in this study were correlated with the observed values, as shown in Figure 11 , and the mean values of carbon storage of both forest types were within the range of the values reported in previous studies ( Table 5 ).

This study shows the spatial differences between drought resistance and resilience of the EBF and ENF in the Chinese subtropics ( Figures 8 , 9 ). As shown in the figure, the resistance value of the EBF to drought was significantly higher than that of the ENF, indicating that EBF has a stronger drought resistance ability in the subtropical region of China, which is consistent with previous studies (Huang and Xia, 2019; Shao et al., 2022). The reason may belong to the higher photosynthesis efficiency of EBF than that of ENF during drought (Wu and Wang, 2022). Previous studies found that EBF could accelerate the loss of old leaves and maintain the growth of young leaves to maintain the light use efficiency, and increase the carbon sequestration capability (Wu et al., 2016; Wei et al., 2017). Although there is significant difference in the resistance of the two forest types, no significant difference in resilience between the two forests 1-4 years after the drought. Previous studies found ecosystem resilience at a large scale can be expressed by their respective WUE levels (Ponce-Campos et al., 2013; Stan et al., 2021). Therefore, due to the similar change trend of WUE in the case of changes in hydroclimatic conditions of EBF and ENF (Sharma and Goyal, 2018; Huang and Xia, 2019), resulting the similar resilience of EBF and ENF to drought in subtropical China. However, the varies in climate and geological conditions cased spatial difference of WUE in two forests, which may be the reason of spatial heterogeneity in resilience throughout subtropical China (Guo et al., 2019; Shao et al., 2022).

EBF and ENF are the predominant vegetation types in subtropical regions of China, especially EBF, which accounts for approximately 60% of global photosynthetic carbon uptake (Mitchard, 2018). Their stability to drought plays a key role in maintaining the stability of subtropical forest ecosystems. According to the analysis results, there were certain differences in the ecosystem stability of the two types of forests, indicating that the stability of different biomes is different (Pennington and Lavin, 2016; Gazol et al., 2018; Anderegg et al., 2020). Current forest management and protection strategies (e.g. Natural forest resources protection, Returning farmland to forest, etc.) have made great improvement on carbon sequestration capability (Kong et al., 2020; Tong et al., 2020), and the increasing forest stand quality further enhanced the LUE and WUE of forest ecosystems (Ponce-Campos et al., 2013; Guo et al., 2019; Huang and Xia, 2019; Shao et al., 2022), leading to strong resistance and resilience to drought. There is no doubt that these strategies should continue to be implemented and widely promoted to provide strong support for subtropical forests to respond and adapt to climate change. However, whether the stability of evergreen forests could remain at current stage under the continuous increasing severe drought events should be further investigated in the future (Easterling et al., 2000; Dai, 2013; Reichstein et al., 2013).