### Investigation of the biogeophysical effect of the idealized deforestation

To obtain the biogeophysical effect of the idealized deforestation, the piControl simulation of the CMIP6^{44} and the idealized global deforestation simulation (deforest-globe) of the LUMIP^{44} are used.

piControl is a long-term (≥500 years) fully coupled simulation and serves as the baseline experiment of CMIP6. piControl is representative of the period prior to the onset of large-scale industrialization, with the reference year being 1850. Volcanic aerosols are fixed at, as closely as possible, the mean values over the historical period of 1850–2014. Solar variability is fixed at the mean value over the first two solar cycles of the historical period (e.g., 1850–1873). Anthropogenic forcings, e.g., GHGs, aerosols, and LULCCs, are fixed at the levels of the year 1850.

deforest-globe is branched from piControl and is an 80-year fully coupled simulation. In deforest-globe, 20 million km^{2} of tree cover is linearly converted to natural grassland cover during the first 50 years. Deforestation is only performed on the top 30% grid cells with the highest initial fractional tree cover (determined by piControl). In the following 30 years, the tree cover is maintained at a constant value to allow the system to equilibrate. All external forcings except the land cover are the same as in piControl. While the response of terrestrial CO_{2} flux to deforestation is allowed, the atmospheric CO_{2} concentration is fixed at the value of piControl. In other words, deforestation-induced CO_{2} emissions do not reach the atmosphere. Thus, deforestation only influences climate through biogeophysical processes.

The last 30-year (year 51–80) simulations since the branching year are used for analysis for piControl and deforest-globe. The biogeophysical effect of the idealized deforestation is obtained by comparing the 30-year simulations of deforest-globe and piControl (deforest-globe minus piControl). By the time of this study, five earth system models contribute to both the piControl and deforest-globe simulations and provide the required daily 2-meter temperature outputs; thus, these five models are used (Supplementary Table S1). For each model, only one member is available for deforest-globe; thus, the single member of deforest-globe and the corresponding member of piControl are used for analysis.

The total tree cover change is obtained by subtracting the tree cover fraction of piControl from the value at year 51 of deforest-globe. Following the same deforestation protocol of deforest-globe, the models are required to produce a spatial deforestation signal that, as closely as possible, replicates that shown in Fig. 1a. Nevertheless, the models still slightly differ in the tree cover changes at regional scales. This intermodel consistency is inevitable and can be attributed to the different initial tree cover fractions and divergent model structures, e.g., vegetation dynamics^{45}. For the models incorporating the dynamic vegetation module, vegetation dynamics are required to be switched off (if possible) on the deforested grid cells and allowed outside the deforested grid cells^{45}. The tree cover of EC-Earth3-Veg and UKESM1-0-LL still evolve during the last 30-year simulation (Fig. 1b) because the vegetation dynamics are switched on outside the deforested grid cells for these two models. Nevertheless, the tree cover change is limited during the last 30-year simulation for these two models.

### Investigation of the biogeophysical effect of historical deforestation

To obtain the biogeophysical effect of deforestation during the historical period, the historical simulation (historical) of CMIP6^{44} and the historical no land-use simulation (hist-noLu) of LUMIP^{45} are used.

historical is a fully coupled simulation covering the historical period of 1850–2014. This simulation is branched from piControl and is driven by evolving and externally imposed forcings. All natural (e.g., volcanic aerosols and solar variability) and anthropogenic (e.g., GHGs, aerosols, and LULCCs) forcings are included and are largely based on observations. The historical LULCCs are prescribed by the second generation of the Land-Use Harmonization dataset^{63} and include detailed changes in land cover, land use, and land management.

hist-noLu is the same as historical except with LULCCs fixed at the levels of the year 1850. Note that the atmospheric CO_{2} concentration in hist-noLu is identical to the value in historical. Thus, the biogeophysical effect of historical deforestation is isolated by comparing the last 30-year (1985–2014) simulations of historical and hist-noLu (historical minus hist-noLu). By the time of this study, nine earth system models contribute to both the historical and hist-noLu simulations and provide the required daily 2-meter temperature output; thus, these nine models are used here (Supplementary Table S1). For each model, multiple members are available for historical and hist-noLu. To save computational cost, the first members of historical and hist-noLu are used for analysis.

The total tree cover change during the historical period is obtained by subtracting the last 30-year (1985–2014) mean tree cover fraction of hist-noLu from the corresponding value of historical. Note that while the historical LULCCs are prescribed by the Land-Use Harmonization dataset, the implementation of the historical LULCCs across different models is commonly different (Supplementary Fig. 7), and this problem has always plagued model intercomparison^{45}. Such intermodel inconsistency may cause the models to diverge in the response of daily temperature variability to the historical LULCCs. Moreover, both tree losses and gains may occur in a region, confounding the signals in daily temperature variability arising from forest losses and gains. However, we identify a large area of forest losses during the historical period in North America (Fig. 5a). This tree cover change is considerable, spatially coherent, and of high intermodel consistency (Supplementary Fig. 7). Therefore, we take North America as a case study to investigate the response of daily temperature variability to historical deforestation.

### Investigation of the net effect of anthropogenic forcings

To obtain the net effect of anthropogenic forcings during the historical period of 1850–2014, the hist-nat simulation of the Detection and Attribution Model Intercomparison Project^{64}, as well as the historical, are used.

hist-nat resembles historical but is forced with only solar and volcanic forcings from historical. Thus, the net effect of all anthropogenic forcings (GHGs, aerosols, and LULCCs) is isolated by comparing the last 30-year (1985–2014) simulations of historical and hist-nat (historical minus hist-nat). By the time of this study, four earth system models contribute to both the historical and hist-nat simulations and provide the required daily 2-meter temperature output; thus, these four models are used here (Supplementary Table S1). For each model, multiple members are available for historical and hist-nat. Similarly, the first members of historical and hist-nat are used for analysis.

### Investigation of the biogeophysical effect of potential afforestation

To obtain the biogeophysical effect of potential afforestation by the end of the twenty-first century, the ssp370 simulation of the Scenario Model Intercomparison Project^{55} and the ssp370-ssp126Lu simulation of the LUMIP are used.

ssp370 is a fully coupled simulation covering the future period of 2015–2100 under the SSP3-7.0 scenario. This scenario represents the medium to high end of the range of future forcing pathways with substantial LULCCs (in particular global forest losses) and high emissions of near-term climate forcers (namely, tropospheric aerosols, tropospheric O_{3} precursors, and CH_{4}). ssp370-ssp126Lu is the same as ssp370, but LULCCs are from the SSP1-2.6 scenario. The SSP1-2.6 scenario represents the low end of the range of future forcing pathways that aim to constrain the global mean temperature below 2 °C by 2100. In contrast to the SSP3-7.0 scenario, there are substantial forest gains due to afforestation or reforestation for the SSP1-2.6 scenario. The LULCCs for both SSP1-2.6 and SSP3-7.0 are prescribed by the Land-Use Harmonization dataset. Note that the atmospheric CO_{2} concentration in ssp370 is identical to the value in ssp370-ssp126Lu. Thus, the biogeophysical effect of potential afforestation is isolated by comparing the last 30-year (2071–2100) simulations of ssp370 and ssp370-ssp126Lu (ssp370-ssp126Lu minus ssp370).

By the time of this study, six earth system models contribute to both the ssp370 and ssp370-ssp126Lu simulations and provide the required daily 2-meter temperature output; thus, these six models are used here (Supplementary Table S1). For each model, only one member is available for ssp370-ssp126Lu. Thus, the single member of ssp370-ssp126Lu and the corresponding member of ssp370 are used for analysis.

### Definitions of the daily temperature variability

To quantify the daily temperature variability, two indices are used: the DTDT^{45} and the SDT^{32}. Specifically, DTDT is defined as the mean value of absolute differences in temperature between every two neighboring days during a given time period, expressed as:

$${{{{{\rm{DTDT}}}}}}=\frac{1}{n-1}\mathop{\sum }\limits_{i=1}^{n-1}\left|{T}_{i+1}-{T}_{i}\right|$$

(1)

where *T*_{i} denotes the 2-meter temperature on day *i*, and *n* denotes the total days of the period.

SDT is a more commonly used index to quantify daily temperature variability and is expressed as:

$${{{{{{\rm{SDT}}}}}}=\sqrt{\frac{1}{n-1}{\sum }_{i=1}^{n}\left({T}_{i}^{*}-\bar{T}\right)^2}}$$

(2)

where \({T}_{i}^{*}\) denotes the 2-meter temperature on day *i* with the annual cycle removed.

Note that DTDT and SDT differ in some aspects, although they both describe daily temperature variability^{31,37,65}. First, SDT is insensitive to the positions of daily temperature anomalies. For example, an orderly (e.g., 25, 25, 25, 25, 15, 15, 15 and 15 °C) and an oscillatory (e.g., 25, 15, 25, 15, 25, 15, 25 and 15 °C) time series of daily temperature have identical SDT values (5.27 °C). However, the temperature changes between neighboring days can be felt once and seven times in the orderly and oscillatory series, respectively. Such a difference between the orderly and oscillatory series can be reflected by DTDT (1.43 versus 10 °C). Second, compared to DTDT, SDT contains variabilities at multiple time scales and may be susceptible to low-frequency temperature variabilities (e.g., weekly variability).

Here, we do not determine which index is better for describing daily temperature variability. Thus, both DTDT and SDT are used to quantify daily temperature variability. More importantly, the results are independent of the metrics of daily temperature variability (Fig. 1c–f versus Supplementary Fig. 4).

Some models provide the daily 2-meter mean temperature that can be used to directly calculate DTDT and SDT. A few models, however, only provide the daily 2-meter maximum and minimum temperatures. For these models, the maximum and minimum temperatures are averaged to produce the estimate of the mean temperature, and then the estimate is further used to calculate DTDT and SDT. The model output variables required to calculate DTDT and SDT are listed in Supplementary Table S2.

### Attribution of the change in daily temperature variability

According to the thermodynamic energy equation^{54}, the temperature difference (δ*T*) between two neighboring days can be expressed as:

$${{{{{\rm{\delta }}}}}}T=\frac{\partial T}{\partial t}=-{{{{{\bf{V}}}}}}\,{\cdot }\, \nabla T+\left(\frac{{RT}}{{c}_{p}P}-\frac{\partial T}{\partial P}\right)\omega+\frac{1}{{c}_{p}}\frac{\partial Q}{\partial t}$$

(3)

where *T* is the daily 2-meter temperature (°C), **V** is the daily 10-meter vector wind (m·s^{−1}), \(\nabla\) is the horizontal divergence operator, *R* is the gas constant for dry air (287 J·K^{−1}·kg^{−1}), *c*_{p} is the specific heat of dry air at constant pressure (1004 J·K^{−1}·kg^{−1}), *P* is the atmospheric pressure (Pa), *ω* is the daily near-surface vertical velocity (Pa·s^{−1}), and *Q* is the daily diabatic heating (J·kg^{−1}).

Substituting Eq. (3) into Eq. (1), Eq. (1) can be rewritten as:

$${{{{{\rm{DTDT}}}}}} \;=\; \frac{1}{n-1}{\sum }_{i=1}^{n-1}\left|{\delta T}_{i}\right |=\frac{1}{n-1}{\sum }_{i=1}^{n-1}\left|{\frac{\partial T}{\partial t}}_{i}\right|\; \approx \; {{{{{\rm{term}}}}}}1+{{{{{\rm{term}}}}}}2+{{{{{\rm{term}}}}}}3$$

(4)

$${{{{{\rm{term}}}}}}1=\frac{1}{n-1}{\sum }_{i=1}^{n-1}\left|-{{{{{{\bf{V}}}}}}}_{i}{\cdot }\nabla {T}_{i}\right|$$

(5)

$${{{{{\rm{term}}}}}}2=\frac{1}{n-1}{\sum }_{i=1}^{n-1}\left|\left(\frac{R{T}_{i}}{{c}_{p}{P}_{i}}-{\frac{\partial T}{\partial P}}_{i}\right){\omega }_{i}\right|$$

(6)

$${{{{{\rm{term}}}}}}3=\frac{1}{n-1}{\sum }_{i=1}^{n-1}\left|\frac{1}{{c}_{p}}{\frac{\partial Q}{\partial t}}_{i}\right|$$

(7)

Then, the deforestation-induced change in DTDT can be expressed as:

$$\Delta {{{{{\rm{DTDT}}}}}}=\Delta {{{{{\rm{term}}}}}}1+\Delta {{{{{\rm{term}}}}}}2+\Delta {{{{{\rm{term}}}}}}3$$

(8)

where \(\Delta\) denotes the difference between the piControl and deforest-globe simulations (deforest-globe minus piControl) or between the historical and hist-noLu simulations (historical minus hist-noLu). Δterm1 denotes the change in the horizontal near-surface TADV. Δterm2 denotes the change in adiabatic compression and vertical advection. Δterm2 is commonly negligibly small in magnitude and can be reasonably omitted^{65}. Δterm3 denotes the change in the daily variability of diabatic heating.

For Δterm3, diabatic heating is considered to be determined by surface SHF because deforestation directly impacts near-surface temperature mainly through the change in SHF^{5,6}. Strictly speaking, to calculate Δterm3, the sensible heat absorbed by the 2-meter air should be used. However, only the total sensible heat that warms the aboveground air is known. To compromise, we use the total sensible heat to qualitatively estimate the contribution of Δterm3 to the DTDT change, expressed as:

$$\Delta {{{{{\rm{term}}}}}}3=\Delta \left(\frac{1}{n-1}\mathop{\sum}\limits_{i=1}^{n-1}\left | \frac{\partial {{{{{\rm{SHF}}}}}}}{\partial t}_{i}\right | \right)=\Delta \left(\frac{1}{n-1}\mathop{\sum}\limits_{i=1}^{n-1}\left|{{{{{{\rm{SHF}}}}}}}_{i+1}-{{{{{{\rm{SHF}}}}}}}_{i}\right|\right)$$

(9)

where SHF_{i} denotes the surface SHF on day *i*. Obviously, Δterm3 is actually the change in the DTDSHF (see Eq. (1)). It should be emphasized that this approach prevents us from showing Δterm3 in a unit equivalent to ΔDTDT. In other words, Δterm3 (Eq. (9)) can only qualitatively explain the DTDT change, e.g., positive Δterm3 leading to increased DTDT.

In summary, the DTDT response to deforestation can be attributed to the change in the horizontal near-surface TADV and the change in the DTDSHF. Note that the effects of changes in albedo, surface roughness, and evaporation efficiency that may arise from deforestation have already been covered by the TADV and DTDSHF effects. For example, the albedo-driven cooling effect and the evapotranspiration-driven warming effect determine the surface temperature change, which directly influences DTDSHF and indirectly influences the TADV through the modification of the horizontal temperature gradient. The lower surface roughness increases near-surface wind, favoring the higher TADV. Similarly, the atmospheric feedbacks^{66,67} that arise from deforestation to surface temperature (e.g., cloud formation and downwards radiation) are also implicitly included in the TADV and DTDSHF effects.

The model output variables required to calculate TADV and DTDSHF are listed in Supplementary Table S2.

### Validation of the model simulation of daily temperature variability

To validate the model behavior in the simulations of DTDT and SDT, two independent reanalysis datasets are used: the ERA5 reanalysis^{46} from the European Center for Medium Weather Forecasting and the NCEP-DOE AMIP-II reanalysis^{47} from the National Centers for Environmental Prediction and National Center for Atmospheric Research. These two reanalysis datasets can provide reliable observations of daily temperature variability^{36} and are widely used in previous studies^{36,38}. The ERA5 reanalysis provides the global daily mean 2-meter temperature at a horizontal resolution of 0.25° (latitude) × 0.25° (longitude) since 1950. The NCEP-DOE AMI-II reanalysis provides global daily mean 2-meter temperature on the T62 Gaussian grid (192 longitude grids × 94 latitude grids) since 1979.

The 30-year (1985–2014) mean DTDT and SDT from the ERA5 reanalysis (Supplementary Figs. 1i–l and 3i–l) and the NCEP-DOE AMIP-II reanalysis (Supplementary Figs. 1m–p and 3m–p) are used as the observational benchmarks. The 30-year (year 51–80) mean DTDT and SDT simulated by piControl (Supplementary Figs. 1a–d and 3a–d) and the 30-year mean (1985–2014) DTDT and SDT simulated by historical (Supplementary Figs. 1e–h and 3e–h) are validated against the reanalysis. Note that the anthropogenic forcings for piControl are fixed at the levels of the year 1850, which are inconsistent with the reanalysis levels. Thus, there might be larger biases in DTDT and SDT simulated by piControl than historical.

Validated with the reanalysis, both piControl and historical can reasonably simulate the global patterns of DTDT and SDT for each season. Statistically speaking, compared to the ERA5 reanalysis, the root mean squared error in the simulated DTDT is between 0.17 and 0.35 °C for piControl and between 0.16 and 0.24 °C for historical across the seasons. When compared to the NCEP-DOE AMIP-II reanalysis, the root mean squared error in DTDT is slightly larger for both piControl (0.25–0.37 °C across the seasons) and historical (0.23–0.47 °C across the seasons). In general, the root mean squared error in the simulated DTDT is highest in winter, followed by spring and autumn, and lowest in summer. A similar result can be obtained for the simulated SDT.

### Validation of the model representation of the biogeophysical effect of deforestation on mean surface temperature

To evaluate the model behavior in the representation of the biogeophysical effect of deforestation on mean surface temperature, a dataset^{48} mapping the biophysical effect of vegetation cover change at the global scale is used. This dataset is established based on the “space-for-time” substitution method and multiple satellite observations. We use the dataset mapping the potential changes in daytime (~13:30 local time) and nighttime (~01:30 local time) surface temperatures as a result of a complete transition from forests to grasslands and croplands. The dataset is provided at monthly time steps and a horizontal resolution of 1° × 1°. The monthly values are averaged over the period of 2008–2012. The daytime and nighttime values are averaged to obtain the estimate of the daily mean value.

Compared to the satellite observations, the models simulate more widespread cooling effects in the northern mid-latitudes (Supplementary Fig. 2). It should be emphasized that such discrepancies cannot be fully attributed to the poor performance of the models. Instead, this discrepancy largely results from the unfair comparison between the observed and simulated biogeophysical effects^{50}. Specifically, the observed biogeophysical effect of deforestation is retrieved based on the “space-for-time” substitution method through the comparison of the spatially adjacent forest and openland pixels. Since the atmosphere is assumed to be largely identical over forest pixels and neighboring openland pixels, the deforestation-induced atmospheric feedbacks to surface temperature (e.g., cloud formation and downwards radiation)^{66,67} are not fully considered. Thus, the observed biogeophysical effect should be considered largely representative of the local effect of deforestation^{15,48}. In contrast, the atmospheric feedbacks to deforestation are fully considered in the coupled models. Thus, the simulated biogeophysical effect includes both the local and nonlocal (through atmospheric feedback) effects of deforestation^{49}. It has been explicitly demonstrated that deforestation causes widespread nonlocal cooling effects in the northern mid- to high-latitudes, with magnitudes comparable to or even larger than the local effects^{49,68}. That is why the simulations indicate a more widespread cooling effect of deforestation in the northern mid-latitudes than the satellite observations^{49,50} (Supplementary Fig. 2). Moreover, the biogeophysical effect of deforestation also depends on the background climate^{69,70}. Thus, the cold bias in the simulated biogeophysical effect can also be partly attributed to the colder background climate of piControl, with the anthropogenic forcings fixed at the levels of the year 1850.

### The paired forest and openland sites

To obtain the observational evidence on the biogeophysical effect of deforestation on daily temperature variability, the observations from 33 flux sites of the FLUXNET2015 Tier 1 dataset^{51} and 7 flux sites of the AmeriFlux dataset^{52} are used. These sites are located in North America and Europe, where the deforestation effect on daily temperature variability is prominent, as suggested by the simulations. These sites provide gap-filled observations of daily 2-meter air temperature, 10-meter wind speed, and surface SHF, which are required in this study. The SHF is observed through the eddy-covariance method and has been corrected to address the surface energy imbalance issue.

The land cover types of these 34 sites include forests, grasslands, croplands, and open shrublands. Following the “space-for-time” substitution logic, a forest site is paired with a neighboring openland (grasslands, croplands, and open shrublands) site. To minimize the influence of the elevation difference between the paired forest and openland sites on our analysis, the pairs with elevation differences exceeding 500 m are discarded. In total, 24 pairs are available for further analysis (Supplementary Fig. 5 and Supplementary Table S3). The mean horizontal distance between the paired sites is 24.6 km, and the mean elevation difference is 71.7 m. The paired sites share almost the same background atmospheric conditions due to the short horizontal distance^{53}. Therefore, any differences in the near-surface atmosphere (e.g., daily temperature variability) between the paired sites can be mostly attributed to the contrasting land cover types. Again, since the “space-for-time” substitution method is used, only the local effects are considered here, whereas the nonlocal effects are mostly not included. For example, the atmosphere above the paired sites tends to be mixed due to horizontal advection. As a result, the effect related to the change in horizontal near-surface TADV is weakened in this case.

### The Kolmogorov–Smirnov two-sample test

The Kolmogorov–Smirnov two-sample test is a nonparametric test to determine if two samples of data are from the same continuous distribution. The null hypothesis is that the two dataset values are from the same continuous distribution. The alternative hypothesis is that these two datasets are from different continuous distributions. The hypothesis test is carried out at the confidence level of 95%. This test is performed using MATLAB.

Support Lumiserver & Cynesys on Tipeee

**Visit
our sponsors**

Wise (formerly TransferWise) is the cheaper, easier way to send money abroad. It helps people move money quickly and easily between bank accounts in different countries. Convert 60+ currencies with ridiculously low fees - on average 7x cheaper than a bank. No hidden fees, no markup on the exchange rate, ever.

**
Now you
can get a free first transfer up to 500£
with your ESNcard. You can access this offer here.**

Source link