# Quantifying the biophysical effects of forests on local air temperature using a novel three-layered land surface energy balance model

Yongxian Su, Liyang Liu, Jianping Wu, Xiuzhi Chen, Jiali Shang, Philippe Ciais, Guoyi Zhou, Raffaele Lafortezza, Yingping Wang, Wenping Yuan, et al.

► To cite this version:

Yongxian Su, Liyang Liu, Jianping Wu, Xiuzhi Chen, Jiali Shang, et al. Quantifying the biophysical effects of forests on local air temperature using a novel three-layered land surface energy balance model. Environment International, 2019, 132, pp.105080. 10.1016/j.envint.2019.105080. hal-02899642

HAL Id: hal-02899642

https://hal.science/hal-02899642

Submitted on 28 Aug 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Quantifying the biophysical effects of forests on local air temperature using
a novel three-layered land surface energy balance model

<figure><img src="image_2.png" /></figure>

Check for
updates

Yongxian Su<sup>a,e</sup>, Liyang Liu<sup>a</sup>, Jianping Wu<sup>a</sup>, Xiuzhi Chen<sup>b,c,*</sup>, Jiali Shang<sup>d</sup>, Philippe Ciais<sup>e</sup>, Guoyi Zhou<sup>f,**</sup>, Raffaele Lafortezza<sup>g,h</sup>, Yingping Wang<sup>i</sup>, Wenping Yuan<sup>b,c</sup>, Yilong Wang<sup>e</sup>, Hongou Zhang<sup>a</sup>, Guangqing Huang<sup>a</sup>, Ningsheng Huang<sup>j</sup>

ARTICLE INFO

Handling Editor: Adrian Covaci

Keywords:

Forest cooling/warming

Forest biophysical effects

Radiative transfer equation

Land-atmosphere model

Local environmental health

Climate warming mitigation

A B S T R A C T

The well-documented energy balance dynamics within forest ecosystems are poorly implemented in studies of the biophysical effects of forests. This results in limitations to the accurate quantification of forest cooling/ warming on local air temperature. Taking into consideration the forest air space, this study proposes a three-layered (canopy, forest air space and soil [CAS]) land surface energy balance model to simulate air temperature within forest spaces (T<sub>af</sub>) and subsequently to evaluate its biophysical effects on forest cooling/warming, i.e., the air temperature gradient (ΔT<sub>a</sub>) between the T<sub>af</sub> and air temperature of open spaces (T<sub>ao</sub>) (ΔT<sub>a</sub> = T<sub>af</sub> - T<sub>ao</sub>). We test the model using field data for 23 sites across 10 cities worldwide; the model shows satisfactory performance with the test data. High-latitude forests show greater seasonal dynamics of ΔT<sub>a</sub>, generating considerable cooling of local air temperatures in warm seasons but minimal cooling or even warming effects during cool seasons, while low-latitude tropical forests always exert cooling effects with less interannual variability. The interannual dynamics of ΔT<sub>a</sub> are significantly related to the seasonality of solar geometry and canopy leaf phenology. The differences between forest canopy temperature (T<sub>c</sub>) and T<sub>ao</sub>, which are the two most important terms attributed by the CAS model in impacting T<sub>af</sub>, explain a large part of forest cooling and warming (May–July: R<sup>2</sup> = 0.35: November–January: R<sup>2</sup> = 0.51). The novel CAS model provides a feasible way to represent the energy balance within forest ecosystems and to assess its impacts on local air temperatures globally.

1. Introduction

Global forest ecosystems can significantly influence land surface
temperatures (Cao et al., 2010; Chen et al., 2012) through their effects
on surface energy balance, i.e., so-called biophysical effects (Lee et al.,
2011; Zeng et al., 2017). Previous satellite- and model-based studies

have had a main focus on large-scale air temperature (Shashua-Bar and Hoffman, 2002; Chang et al., 2007; Shashua-Bar et al., 2009; Cao et al., 2010; Lee et al., 2011; Chen et al., 2012; Feyisa et al., 2014) or soil surface temperature (Bonan, 2008; Peng et al., 2014; Zhang et al., 2014; Li et al., 2015; Alkama and Cescatti, 2016), and very few have investigated the local air temperature within forest ecosystems.

<footnote><sup>a</sup>Key Lab of Guangdong for Utilization of Remote Sensing and Geographical Information System, Guangdong Open Laboratory of Geospatial Information Technology and Application, Guangzhou Institute of Geography, Guangzhou 510070, China</footnote>

<footnote><sup>b</sup> Guangdong Province Key Laboratory for Climate Change and Natural Disaster Studies, School of Atmospheric Sciences, Sun Yat-sen University, Guangzhou 510275, China</footnote>

<footnote><sup>c</sup> Southern Marine Science and Engineering Guangdong Laboratory, Zhuhai 519082, China</footnote>

<footnote><sup>d</sup> Ottawa Research and Development Centre, Agriculture and Agri-Food Canada, Ottawa, Ontario, Canada</footnote>

<footnote><sup>e</sup> Laboratoire des Sciences du Climat et de l'Environnement, UMR 1572 CEA-CNRS UVSQ, 91191 Gif sur Yvette, France</footnote>

<footnote><sup>f</sup> Institute of Ecology, Jiangsu Key Laboratory of Agricultural Meteorology, Nanjing University of Information Science & Technology, Nanjing 210044, China</footnote>

<footnote><sup>g</sup> Department of Geography, The University of Hong Kong, Centennial Campus, Pokfulam Road, Hong Kong</footnote>

<footnote><sup>h</sup> Department of Agricultural and Environmental Sciences, University of Bari “A. Moro”, Via Amendola 165/A, 70126 Bari, Italy</footnote>

<footnote><sup>i</sup> CSIRO Oceans and Atmosphere, Aspendale, Victoria 3195, Australia</footnote>

<footnote><sup>j</sup> Guangzhou Institute of Geochemistry, Chinese Academy of Sciences, Guangzhou 510640, China</footnote>

<footnote>* Correspondence to: X. Chen, Guangdong Province Key Laboratory for Climate Change and Natural Disaster Studies, School of Atmospheric Sciences, Sun Yat-sen University, Guangzhou 510275, China.</footnote>

<footnote>** Corresponding author.</footnote>

<footnote>E-mail addresses: chenxzh73@mail.sysu.edu.cn (X. Chen), gyzhou@nuist.edu.cn (G. Zhou).</footnote>
Acronyms

$T_{af}$ air temperature of forest space

$T_{ao}$ air temperature of open space

$\Delta T_a$ gradient between $T_{af}$ and $T_{ao}$

$T_s$ soil surface temperature

$T_c$ canopy surface temperature

$\emptyset_n$ downward solar shortwave radiation

$LW$ surface longwave radiation

$LW_{sky}$ downward sky longwave radiation

$LW_{canopy}$ canopy longwave radiation

$LW_{canopy, \downarrow}$ downward canopy longwave radiation

$LW_{soil}$ upward soil longwave radiation

$LE$ total surface latent heat flux

$LE_{soil}$ latent heat loss of soil layer

$LE_{canopy}$ latent heat loss of canopy layer

$H$ total surface sensible heat flux

$H_{soil \to af}$ sensible heat loss of soil layer

$H_{af \to canopy}$ sensible heat flux between forest air space and tree canopy

$H_{canopy \to ao}$ sensible heat flux between canopy and open-air space

$G_{soil}$ heat storage in soil

$G_{tree}$ heat storage in tree

$ET$ evapotranspiration

$c$ extinction coefficient

$LAI$ leaf area index

$u$ cosine of the zenith angle ($\theta$) of the Sun

$\Gamma$ canopy transmission of direct radiation

$K$ von Karman's constant

$\delta = \frac{1 + K\mu}{1 - K\mu}$

$\sigma$ Stefan-Boltzmann constant

$\rho_a$ density of air

$C_p$ specific heat capacity of air

$\varepsilon_s$ proportion of soil surface emissivity

$\varepsilon_{sky}$ proportion of sky emissivity

$C_{cover}$

$r$ aerodynamic resistance to sensible heat transfer

$r_s$ $r$ between soil and forest air spaces

$r_{a, c}$ $r$ between forest air spaces and canopy

$r_{c, a}$ $r$ between canopy and open-air spaces

$\partial$ derivative of $T_{af}$ to environment variables

$C$ relative contributions of environment variables to $T_{af}$

$VPD$ vapor pressure deficit

$h_c$ canopy height

$U(V_z)$ wind speed at reference height

Currently, there has been growing interest in the impact of land cover changes on human health via increased heat exposure (Wolff et al., 2018). Given that the forest canopy covers > 30% of global land surface (FAO, 2016), air temperature within the forest space ($T_{af}$) plays an important role in mitigating local climate warming (Oke, 1989; Bounoua et al., 2002; Georgi and Zafiriadis, 2006; Chang et al., 2007; Shashua-Bar et al., 2010; Oliveira et al., 2011; Susca et al., 2011; Chen et al., 2012); it is all the more important for maintaining local environmental health (Wolff et al., 2018) as well as ecosystem stability (Ellison et al., 2017).

Forest cooling and warming on the local environment, i.e., the gradients between $T_{af}$ and air temperature of open space ($T_{ao}$) ($\Delta T_a = T_{af} - T_{ao}$), have been extensively studied using time-series field observations (Georgi and Zafiriadis, 2006; Potchter et al., 2006; Chang et al., 2007; Zhang et al., 2013). However, previous studies used only a few sites within a single region (Wong and Yu, 2005; Bowler et al., 2010; Feyisa et al., 2014); therefore, their conclusions are likely regional, and may not be applicable elsewhere (Marland et al., 2003; Anderson et al., 2011; Li et al., 2015; Zeng et al., 2017). For an

example, the directions and magnitudes of $\Delta T_a$ reported in the literature differ greatly, ranging from -5.6 °C to +3.3 °C (Jauregui, 1991; Taha et al., 1991; Wong and Yu, 2005; Chang et al., 2007; Potchter et al., 2012). Additionally, multiple factors, e.g., evapotranspiration (ET) (Li et al., 2015; Devaraju et al., 2018), aerodynamic roughness (Lee et al., 2011; Burakowski et al., 2018; Devaraju et al., 2018), albedo (Lee et al., 2011; Devaraju et al., 2018), surface resistance (Li et al., 2019) or background climate (Pitman et al., 2011), were identified as dominant contributors in the satellite- and model-based studies that address large-scale air temperature and soil surface temperature. Whether this mechanism still stands in driving the cooling and warming effects of forests on local air temperatures deserves further in-depth study.

The energy balance dynamics within forest ecosystems have been well documented in previous studies (Choudhury and Monteith, 1988; Monteith and Unsworth, 2007). However, these processes are poorly implemented in relative studies of modeling biophysical effects of forests on local air temperature. Taking into consideration the air space inside a forest ecosystem (hereafter referred to as the forest air space), we divided the forest land into three vertical layers: canopy, forest air

<figure><img src="image_3.png" /><figcaption>Fig. 1. Illustration of the three-layered CAS three-layered land surface energy balance model.</figcaption></figure>
space, and soil surface (CAS) (Fig. 1). The objectives of this study were to: (1) develop a three-layered CAS land surface energy balance model to simulate $T_{af}$; (2) use the CAS model to quantify forest cooling and warming effects (i.e., using $T_{af}$ minus $T_{ao}$ to estimate $\Delta T_a$) on the local environment; and (3) use the CAS model to attribute relative contributions of major climatic drivers to $T_{af}$ based on variable sensitivity analysis. Both the modeled $T_{af}$ and derived $\Delta T_a$ were validated using a global dataset of 373 field observations collected at 23 sites in 10 cities (Fig. 2, Supplementary dataset).

## 2. Methods, model development and datasets

### 2.1. The novel three-layered CAS land surface energy balance model

#### 2.1.1. Model development strategies

The classic land surface energy balance model (Norman et al., 1995; Ca et al., 1998; Monteith and Unsworth, 2007; Wu et al., 2007; Jones and Rotenberg, 2011; Vidrih and Medved, 2013) treats the forest land as one composite layer and does not isolate the land surface layer into multiple layers (e.g., the canopy, forest air space and soil layers).

In this study, we divided the forest land surface into three vertical layers: canopy, forest air space, and soil surface (left side in Fig. 1). In contrast to the classic land surface energy balance model (Supplementary methods), which treats the forest land surface as an intact entity, the land surface energy balances in this study are partitioned into soil surface energy balance and canopy energy balance (right side in Fig. 1). We initially developed the soil surface energy balance (Section 2.1.2) and canopy energy balance (Section 2.1.3) equations, followed by the three-layered CAS land surface energy balance equation (Section 2.1.4). Lastly, we derived the $T_{af}$ estimation model (Section 2.1.5) from the novel CAS model.

#### 2.1.2. Soil surface energy balance

For a soil surface underneath the canopy, the incoming heat fluxes (left part of Eq. (1)) are mainly from three sources: the downward solar shortwave radiation ($\varnothing_n$) and downward sky longwave radiation ($LW_{sky}$) penetrating through the tree canopy, and the downward canopy longwave radiation ($LW_{canopy, \downarrow}$). We considered the transmission ($\Gamma$, the exponential term in Eq. (1)) as the penetration rate of $\varnothing_n$ and $LW_{sky}$ that reached the understory of the soil surface. Similarly, $LW_{canopy, \downarrow}$ is defined as the product of canopy longwave radiation ($LW_{canopy}$) and $\Gamma$,

i.e., $LW_{canopy, \downarrow} = \Gamma L W_{canopy}$. $\Gamma$ is the transmission of the direct radiation through the canopy layers (He et al., 2017).

The outflow of energy fluxes from the understory soil surface occurs mainly in three ways (right part of Eq. (1)): the latent heat loss of soil layer ($LE_{soil}$), the sensible heat loss of soil layer ($H_{soil \to af}$), the upward longwave radiation of soil layer ($LW_{soil}$), and the heat storage in soil ($G_{soil}$) (Supplementary methods). The energy balance equation of the understory soil surface is:

$$ (\varnothing_n + LW_{sky} + LW_{canopy}) \exp\left(-\frac{cLAI}{u}\right) = LE_{soil} + LW_{soil} + H_{soil \to af} + G_{soil} \quad (1) $$

where $LAI$ is the leaf area index, $c$ is the extinction coefficient (yielding a value of 0.5), and $u$ is the cosine of the zenith angle ($\theta$) of the Sun (He et al., 2017). The exponential term $\exp(-\frac{cLAI}{u})$ is the transmission $\Gamma$. $cLAI$ is the canopy optical thickness (He et al., 2017).

#### 2.1.3. Canopy energy balance

The energy fluxes absorbed by the tree canopy (left part of Eq. (2)) are $\varnothing_n$, $R_{sky}$, and $R_{soil}$. The term $(1 - \Gamma)$ in the tree canopy part is the interception proportion of radiation ($\varnothing_n$, $LW_{sky}$, and $LW_{soil}$), which includes the absorption and reflection components, given as $\varnothing_n [1 - \exp(-\frac{cLAI}{u})] \cdot LW_{sky} [1 - \exp(-\frac{cLAI}{u})] \text{ and } LW_{soil} [1 - \exp(-\frac{cLAI}{u})]$, respectively.

The energy flux loss from the tree canopy mainly goes through five paths (right part of Eq. (2)): latent heat loss resulting from forest canopy transpiration ($LE_{canopy}$), the sensible heat flux between forest air space and tree canopy ($H_{af \to canopy}$), the sensible heat flux between canopy and open air space ($H_{canopy \to ao}$), $LW_{canopy}$ flux, and the heat storage in tree ($G_{tree}$).

$$ (\varnothing_n + LW_{sky} + LW_{soil}) \left[ 1 - \exp\left(-\frac{cLAI}{u}\right) \right] = LE_{canopy} + H_{af \to canopy} + H_{canopy \to ao} + LW_{soil} + G_{tree} \quad (2) $$

#### 2.1.4. The three-layered CAS land surface energy balance model

By combining the soil surface energy balance equation (Eq. (1)) and canopy energy balance equation (Eq. (2)), we ultimately derived the three-layered CAS land surface energy balance equation (Eq. (3)).

<figure><img src="image_4.png" /><figcaption>Fig. 2. Locations of the 23 sites collected from the published literatures. Northern high latitudes: > 23.5°N; Tropics: 23.5°S-23.5°N; Southern high latitudes: > 23.5°S.</figcaption></figure>
$$ 
\begin{aligned}
\emptyset_n + LW_{\text{sky}} &= LE + H_{\text{soil→af}} + H_{\text{af→canopy}} + H_{\text{canopy→ao}} + \\
& \qquad \left\{ LW_{\text{soil}} \exp\left(-\frac{cLAI}{u}\right) + LW_{\text{canopy}} \left[1 - \exp\left(-\frac{cLAI}{u}\right)\right] \right\} + G_{\text{soil}} + G_{\text{tree}}
\end{aligned}
\quad (3) $$

where LE is the total latent heat flux of the forest land surfaces (LE = LE<sub>canopy</sub> + LE<sub>soil</sub>); G<sub>tree</sub> ≈ 0 (Allen et al., 1998), as this part of the fluxes into the forest is much smaller than the other parts.

By partitioning the forest land surface into canopy, forest air space and soil surface, the CAS land surface energy balance model decomposes the surface sensible heat flux (H) into three parts – H<sub>soil→af</sub>, H<sub>af→canopy</sub> and H<sub>canopy→ao</sub> and decomposes the surface longwave radiation (LW) into two parts – LW<sub>soil</sub> exp(-cLAI/u) and LW<sub>canopy</sub> [1 - exp(-cLAI/u)]. Compared with the classic land surface energy balance equation, the three-layered CAS land surface energy balance model is a more accurate representation of the energy flux dynamics occurring within forest ecosystems. The formulas to compute the above various fluxes are introduced in the Supplementary methods.

## 2.1.5. The T<sub>af</sub> retrieval model deduced from the novel CAS model

The T<sub>af</sub> estimation model (Eq. (4)) was derived by forcing the calculation functions of LW<sub>sky</sub>, LW<sub>soil</sub>, LW<sub>canopy</sub>, H<sub>soil→af</sub>, H<sub>af→canopy</sub>, H<sub>canopy→ao</sub>, and G<sub>soil</sub> into the CAS land surface energy balance equation. For a detailed step-by-step deduction, please refer to the Supplementary methods. As can be seen from Eq. (4), T<sub>af</sub> is impacted jointly by ∅<sub>n</sub>, LE, soil surface temperature (T<sub>s</sub>), canopy surface temperature (T<sub>c</sub>), and T<sub>ao</sub> under different LAI and u.

$$ 
\begin{aligned}
\left(\frac{1}{r_s} - \frac{1}{r_{a,c}}\right) \rho_a C_p T_{af} \\
= -\frac{1}{\delta} \left[ 1 + \exp\left(-\frac{cLAI}{u}\right) \right] \emptyset_n + \left[ \frac{1}{\delta} + \exp\left(-\frac{cLAI}{u}\right) \right] LE + \frac{\rho_a C_p}{\delta} \\
\left(\frac{1}{r_{c,a}} - \frac{1}{r_{a,c}}\right) T_c + \frac{\rho_a C_p}{r_s} T_s + \frac{1}{\delta} \left[ 1 + \exp\left(-\frac{cLAI}{u}\right) \right] \varepsilon_s \sigma T_s^4 - \frac{\rho_a C_p}{r_{c,a}} \\
T_{ao} - \frac{1}{\delta} \left[ 1 + \exp\left(-\frac{cLAI}{u}\right) \right] \varepsilon_{sky} \sigma [\varepsilon_{sky} + 0.8(1 - \varepsilon_{sky}) C_{cover}]^4 T_{ao}^4
\end{aligned}
\quad (4) $$

where r<sub>s</sub>, r<sub>a,c</sub> and r<sub>c,a</sub> (unit: S m<sup>-1</sup>) are the resistance to sensible heat transfer between soil and forest air spaces, forest air spaces and canopy, and canopy and open air spaces, respectively; ρ<sub>a</sub> is the density of air (ρ<sub>a</sub> = 1.25 kgm<sup>-3</sup>); C<sub>p</sub> is the specific heat capacity of air (C<sub>p</sub> = 1004 Jkg<sup>-1</sup>k<sup>-1</sup>); δ = (1 + K<sub>μ</sub>) / (1 - K<sub>μ</sub>), K is von Karman's constant, producing a value of 0.41; ε<sub>s</sub> and ε<sub>sky</sub> are the proportions of soil surface and sky emissivity (ε<sub>s</sub> = 0.81; ε<sub>sky</sub> = 0.9), respectively; σ is the Stefan-Boltzmann constant (5.67 * 10<sup>-8</sup> Wm<sup>-2</sup>K<sup>-4</sup>); and C<sub>cover</sub> is the cloud coverage rate.

<table><caption>**Table 1**
Global gridded datasets of ∅<sub>n</sub>, T<sub>ao</sub>, VPD, ET, T<sub>s</sub>, LAI, U(V<sub>z</sub>), C<sub>cover</sub> and h<sub>c</sub> products.</caption><thead><tr><th>Variable</th><th>Source/name</th><th>Version</th><th>Spatial solution</th><th>Temporal solution</th><th>Reference</th></tr></thead><tbody><tr><td>∅<sub>n</sub></td><td>Climatic Research Unit (CRU) downward shortwave radiation</td><td>-</td><td>0.5*</td><td>6-Hour</td><td>Wei et al. (2014)</td></tr><tr><td>T<sub>ao</sub></td><td>Climatic Research Unit (CRU) open air temperature</td><td>-</td><td>0.5*</td><td>6-Hour</td><td>Wei et al. (2014)</td></tr><tr><td>VPD</td><td>Climatic Research Unit (CRU) vapor pressure deficit</td><td>4.01</td><td>0.5*</td><td>Monthly</td><td>Harris et al. (2017)</td></tr><tr><td>ET</td><td>Global Terrestrial Evapotranspiration Data Set (MOD16A2)</td><td>-</td><td>1 km</td><td>8-Day</td><td>Mu et al. (2011)</td></tr><tr><td>T<sub>s</sub></td><td>MODIS Land Surface Temperature and Emissivity (MYD11C1)</td><td>Collection 6</td><td>0.05*</td><td>Daily</td><td>Wan (2014)</td></tr><tr><td>LAI</td><td>MCD15A3H MODIS/Terra + Aqua Leaf Area Index/FPAR (MOD15A3H)</td><td>Collection 6</td><td>500 m</td><td>4-Day</td><td>Myneni et al. (2002)</td></tr><tr><td>U(V<sub>z</sub>)</td><td>ERA Interim wind speed</td><td>-</td><td>~79 km</td><td>6-Hour</td><td>Berrisford et al. (2011)</td></tr><tr><td>C<sub>cover</sub></td><td>NOAA Climate Cloud fraction Data Record (CDR)</td><td>5.3</td><td>0.1*</td><td>Daily</td><td>Heidinger et al. (2014)</td></tr><tr><td>h<sub>c</sub></td><td>Global 1 km Forest Canopy Height</td><td>-</td><td>1 km</td><td>-</td><td>Simard et al. (2011)</td></tr></tbody></table>

## 2.2. Sensitivity analysis to decompose the contributions of underlying drivers to T<sub>af</sub>

### 2.2.1. Sensitivity analysis of T<sub>af</sub> to influencing factors

The derivative method has been widely used in climate-related sensitivity analysis to test the impact of nonlinear independent variables on dependent variables (Paltridge, 1980; Dooge, 1992; Friedlingstein et al., 2003; Zhou et al., 2015). Herein, we used the derivatives of T<sub>af</sub> to ∅<sub>n</sub>(∂∅<sub>n</sub>), LE(∂LE), T<sub>s</sub>(∂T<sub>s</sub>), T<sub>c</sub>(∂T<sub>c</sub>) and T<sub>ao</sub>(∂T<sub>ao</sub>) to evaluate the corresponding variable sensitivity. The functions of ∂∅<sub>n</sub>, ∂LE, ∂T<sub>s</sub>, ∂T<sub>c</sub>, and ∂T<sub>ao</sub> are given individually as:

$$ \partial_{\emptyset_n} = - \frac{\left[ 1 + \exp\left(-\frac{cLAI}{u}\right) \right]}{\left(\frac{1}{r_s} - \frac{1}{r_{a,c}}\right) \delta \rho_a C_p} \quad (5) $$

$$ \partial_{LE} = \frac{\frac{1}{\delta} + \exp\left(-\frac{cLAI}{u}\right)}{\left(\frac{1}{r_s} - \frac{1}{r_{a,c}}\right) \rho_a C_p} \quad (6) $$

$$ \partial_{T_c} = \frac{\frac{\rho_a C_p}{\delta} \left( \frac{1}{r_{c,a}} - \frac{1}{r_{a,c}} \right)}{\left( \frac{1}{r_s} - \frac{1}{r_{a,c}} \right)} \quad (7) $$

$$ \partial_{T_s} = \frac{\frac{4}{\delta} \varepsilon_s \sigma \left[ 1 + \exp\left(-\frac{cLAI}{u}\right) \right] T_s^3 + \frac{\rho_a C_p}{r_s}}{\left(\frac{1}{r_s} - \frac{1}{r_{a,c}}\right) \rho_a C_p} \quad (8) $$

$$ \partial_{T_{ao}} = -\frac{\frac{4}{\delta} \left[ 1 + \exp\left(-\frac{cLAI}{u}\right) \right] \varepsilon_{sky} \sigma [\varepsilon_{sky} + 0.8(1 - \varepsilon_{sky}) C_{cover}] T_a^3 + \frac{\rho_a C_p}{r_{c,a}}}{\left(\frac{1}{r_s} - \frac{1}{r_{a,c}}\right) \delta \rho_a C_p} \quad (9) $$

### 2.2.2. Relative contributions of the influencing factors to T<sub>af</sub>

Zhou et al. (2015) defined the relative contribution as the rate of individual product of each variable and its sensitivity to the sum of the products of all variables and their corresponding sensitivity. Note that Zhou et al.'s (2015) method was targeted at evaluating general mean contributions of independent variables (X) to a dependent variable (Y) (Zhou et al., 2018), which is completely different from the time-series studies that calculated contributions of X variances to Y variances. The former one decomposes total magnitudes of Y into X variables, while the latter one relates the magnitude changes of Y to changes of X variables.

In this study, we adopted Zhou et al.'s (2015) method to calculate the general mean contributions of ∅<sub>n</sub>(C<sub>∅n</sub>), LE(C<sub>LE</sub>), T<sub>s</sub>(C<sub>Ts</sub>), T<sub>c</sub>(C<sub>Tc</sub>), and T<sub>ao</sub>(C<sub>Tao</sub>) to T<sub>af</sub>. The sum of the products of the ∅<sub>n</sub>, LE, T<sub>s</sub>, T<sub>c</sub>, and T<sub>ao</sub> variable terms and their corresponding sensitivities are given in Eq. (10). The formulas of C<sub>∅n</sub>, C<sub>LE</sub>, C<sub>Ts</sub>, C<sub>Tc</sub>, and C<sub>Tao</sub> are expressed in Eqs. (11)-(15), respectively.
$$SUM = |\partial_{\emptyset_n} \emptyset_n| + |\partial_{LE} LE| + |\partial_{T_c} T_c| + |\partial_{T_s} T_s| + |\partial_{T_{ao}} T_{ao}| \quad (10)$$

$$C_{\emptyset_n} = \left| \frac{\partial_{\emptyset_n} \emptyset_n}{SUM} \right| * 100\% \quad (11)$$

$$C_{LE} = \left| \frac{\partial_{LE} LE}{SUM} \right| * 100\% \quad (12)$$

$$C_{T_c} = \left| \frac{\partial_{T_c} T_c}{SUM} \right| * 100\% \quad (13)$$

$$C_{T_s} = \left| \frac{\partial_{T_s} T_s}{SUM} \right| * 100\% \quad (14)$$

$$C_{T_{ao}} = \left| \frac{\partial_{T_{ao}} T_{ao}}{SUM} \right| * 100\% \quad (15)$$

## 2.3. Model inputs and validations

### 2.3.1. Model input datasets for global application

Multi-year global gridded datasets – i.e., the National Centers for Environmental Prediction (NCEP) and Climatic Research Unit-NCEP (CRUNCEP) $\phi_n$ products, CRUNCEP $T_{ao}$ products, CRUNCEP Vapor pressure deficit (VPD) products, Moderate Resolution Imaging Spectroradiometer (MODIS) ET products, MODIS daytime $T_s$ products, and MODIS LAI products (Table 1) – were averaged as input data to force Eq. (4) to simulate global seasonal $T_{af}$ and then to estimate $\Delta T_a$ using simulated $T_{af}$ minus CRUNCEP $T_{ao}$. Similarly, the above datasets were also used as inputs to calculate variable sensitivity ($\partial_{\emptyset_n}, \partial_{LE}, \partial_{T_c}, \partial_{T_s}$ and $\partial_{T_{ao}}$, Eqs. (5)–(9)) and relative contributions ($C_{\emptyset_n}, C_{LE}, C_{T_c}, C_{T_s}$ and $C_{T_{ao}}$, Eqs. (11)–(15)), respectively.

To better match the model outputs with the field observations, we extracted input data from the finest spatial and temporal resolutions of raster imageries (Table 1) according to the temporal and location records of field experiments (Supplementary dataset). For global simulations, we resampled the raster imageries of input data into monthly temporal resolution with 0.5° × 0.5° spatial resolution using the bilinear interpolation method (Liu et al., 2017).

### 2.3.2. Field observations for model validation

For model validations, 373 field samples from 23 sites (Fig. 2) in 10 cities across tropical and temperate latitudes were collected from globally published literature (see Supplementary dataset). Of the 23 sites, four are tropical evergreen forests and 19 are located at northern high latitudes (17 deciduous and 2 evergreen forests). These data contain the following information: author information, city, latitude, longitude, year, month, day, and observed $T_{af}$ and $\Delta T_a$.

The supplementary dataset also contains some variables extracted from the global gridded products, such as tree canopy height ($h_c$), wind speed at reference height $U(V_z)$, cloud coverage rate ($C_{cover}$), LAI, daytime $T_s$, ET, $\emptyset_n$ and $T_{ao}$. The information of the raster datasets employed is listed in Table 1. Other variables, i.e., $\mu$, $r_s$, $r_a$, $c$, $r_c$, $a$, and $T_c$, are estimated. Detailed methods for estimating these parameters are described in the Supplementary methods. The cross-validation results between the CAS model retrieved $T_c$ and the forward simulated $T_c$ from $T_a$, $\emptyset_n$, and VPD are shown in Fig. S1. The simulated $T_{af}$ was derived using Eq. (4), and simulated $\Delta T_a$ was calculated using simulated $T_{af}$ minus CRUNCEP $T_{ao}$.

## 3. Results

### 3.1. Evaluations of simulated $T_{af}$ and $\Delta T_a$ using globally published data

In total, 373 field data records containing $T_{af}$ and $\Delta T_a$ (Section 2.3.2, Supplementary dataset) were collected from the published literature to validate the CAS model outputs. Of the 373 samples, the lowest value of field $T_{af}$ (+5.5 °C) was observed in the winter season of Beijing located at a northern high latitude (longitude: 116.42°E, latitude: 40.05°N) (Ji et al., 2012); the highest $T_{af}$ (+42.1 °C) was also observed in Beijing (longitude: 116.39°E, latitude: 39.37°N), but during the hot summer season (Liu et al., 2008). Apparently, high-latitude forests exhibit greater seasonal variability of $T_{af}$ compared with low-latitude tropical forests. Despite the considerable variation that exists in the collected $T_{af}$ (Fig. 3a), the proposed CAS model succeeded in estimating seasonal $T_{af}$ across different climatic zones (tropical and northern/southern high latitudes), with a root mean square error (RMSE) of 0.98 °C (R² = 0.98) (Fig. 3c). Accounting for the potential impacts of background temperatures (e.g., heat island effect) on the CAS model, we plotted the scatter diagrams between the $T_{ao}$ and $T_{af}$ estimation errors (Fig. S2).

<figure><img src="image_5.png" /><figcaption>Fig. 3. Field data records of $T_{af}$ and $\Delta T_a$ from the published literature (N = 373, Supplementary dataset) for model evaluations. (a–b) Distributions of the collected $T_{af}$ and $\Delta T_a$; (c–d) Scatter plots between simulated $T_{af}$, $\Delta T_a$ and field observations.</figcaption></figure>
Results show that there was no significant trend in T<sub>af</sub> estimation errors with T<sub>ao</sub> changes (R<sup>2</sup> = 0.03), suggesting no systematic impact on the simulation accuracy of the CAS model from background climate. The observed ΔT<sub>a</sub> (Fig. 3b) ranged from -5.6 °C to +2.1 °C. Approximately 87.4% of the measurements revealed a cooling effect, 9.7% a warming

effect, and the remaining 2.9% showed no significant effect. A comparison of these results with field observations showed that the total RMSE of the simulated ΔT<sub>a</sub> was 0.98 °C (R<sup>2</sup> = 0.36) (Fig. 3d) and that 62.5% of absolute errors of the modeled ΔT<sub>a</sub> was lower than 1.0 °C.

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Global patterns of simulated seasonal T<sub>af</sub> and ΔT<sub>a</sub>. (a) Simulated T<sub>af</sub> in the May–July period; (b) Simulated T<sub>af</sub> in the November–January period; (c) Statistics of pixels with different T<sub>af</sub> values at the northern high latitudes (N 23.5°N), tropics (23.5°S–23.5°N) and southern high latitudes (S 23.5°S); (d) Simulated ΔT<sub>a</sub> in the May–July period; (e) Simulated ΔT<sub>a</sub> in the November–January period; (f) Statistics of pixels with different ΔT<sub>a</sub> values at the northern high latitudes (N 23.5°N), tropics (23.5°S–23.5°N) and southern high latitudes (S 23.5°S).</figcaption></figure>
## 3.2. Global patterns of simulated $T_{af}$ and $\Delta T_a$ at different latitudes

The global simulated $T_{af}$ in the two seasonal periods May–July and November–January is shown in Fig. 4a and b. Among the three latitude zones (Fig. 4c and Table 2), the forests at the northern high latitudes had the lowest $T_{af}$ during both the May–July (most frequently occurring $T_{af}$ = 9.13 °C; proportions of $0 < T_{af} < 10$ °C: 58.5%) and November–January (most frequently occurring $T_{af}$ = -2.75 °C; proportions of $T_{af} < 0$ °C: 66.5%) periods, followed by the forests at the southern high latitudes (May–July: most frequently occurring $T_{af}$ = 7.29 °C, proportions of $0 < T_{af} < 10$ °C: 74.0%; November–January: most frequently occurring $T_{af}$ = 21.96 °C, proportions of $T_{af} > 20$ °C: 47.5%). As is typical of tropical forests, $T_{af}$ was the highest in both the May–July (most frequently occurring $T_{af}$ = 20.28 °C; proportions of $10 < T_{af} < 20$ °C: 57.3%) and November–January (most frequently occurring $T_{af}$ = 29.73 °C; proportions of $T_{af} > 20$ °C: 82.3%) periods.

The global estimated seasonal $\Delta T_a$, which was calculated using simulated $T_{af}$ minus $T_{ao}$ (see Fig. 4d and e), ranged from -13.5 °C to +4.9 °C with large seasonal variation and spatial heterogeneity. The pixel number and proportion statistics of $\Delta T_a$ under different thresholds are shown in Fig. 4f and Table 3, respectively. During the May–July period, about 97.9% of forest-induced $\Delta T_a$ was negative at global scale (northern high latitudes: 99.3%; tropics: 97.5%; southern high latitudes: 85.1%), indicating significant net cooling effects of forest ecosystems on the local environment. The most frequently occurring $\Delta T_a$ was estimated to be -4.25 °C at the northern high latitudes, -2.22 °C in the tropics, and -2.59 °C at the southern high latitudes. Results indicate that the forests at northern high latitudes tended to generate larger quantities of negative $\Delta T_a$ than in the low-latitude tropical forests in the May–July seasonal period. During the November–January period (Fig. 4e), however, the forests at the northern high latitudes were more inclined to exhibit warming effects (mostly frequently occurring at $\Delta T_a$ = +1.11 °C), while the dense-canopy forests in low-latitude tropics continued to show cooling effects (mostly frequently occurring at $\Delta T_a$ = -0.20 °C), but at a much reduced magnitude in comparison with the May–July period. On the contrary, forests at the southern high latitudes showed consistent effects as did the tropics, but at a much larger magnitude (most frequently occurring at $\Delta T_a$ = -8.37 °C) due to the November–January warm seasons in the southern hemisphere. In general, the average $\Delta T_a$ at the northern and southern high latitudes exhibited greater seasonal dynamics compared with $\Delta T_a$ at the low tropical latitudes.

## 3.3. Solar geometry and canopy phenology mediate the seasonality of $T_{af}$ and $\Delta T_a$

Theoretically, the three-layered CAS land surface energy balance model depicts five major variables impacting $T_{af}$ (i.e., $\emptyset_n$, LE, $T_s$, $T_c$, and $T_{ao}$), which are jointly mediated by solar geometry (i.e., u) and canopy phenology (i.e., LAI) (Eq. (4)). The sensitivities of $T_{af}$ to $\emptyset_n$, LE, $T_s$, $T_c$ and $T_{ao}$ are represented by the derivatives $\partial_{\emptyset_n}$, $\partial_{LE}$, $\partial_{T_s}$, $\partial_{T_c}$, and $\partial_{T_{ao}}$, respectively (Eqs. (5)–(9)), and plotted in Figs. 5 and 6. In Fig. 5a, $\partial_{\emptyset_n}$ exhibited a decreasing trend as LAI increased. However, the

relationship between $\partial_{\emptyset_n}$ and LAI differed greatly under different $\theta$. This is also reflected in Fig. 6a, where the curve is unimodal when LAI < 1.0 (the peak values of $\partial_{\emptyset_n}$ fell within the range of 70° < $\theta$ < 90°). The $\partial_{T_s}$ exhibited a similar trend along with LAI and $\theta$ as those of $\partial_{\emptyset_n}$, but the $\partial_{T_s}$ values were negative (Figs. 5c and 6c). The $\partial_{LE}$ decreased following a logistic curve as LAI increased (Fig. 5b). $\partial_{T_c}$ and $\partial_{T_{ao}}$ exhibited a slight or no response to the change in LAI (Fig. 5d and e). The $\partial_{LE}$, $\partial_{T_s}$, and $\partial_{T_{ao}}$ all exhibited an exponential relationship with $\theta$ (Fig. 6c, d and e). Overall, $\partial_{\emptyset_n}$ and $\partial_{T_{ao}}$ were always positive, while $\partial_{LE}$, $\partial_{T_s}$, and $\partial_{T_c}$ were always negative. In addition, $\partial_{\emptyset_n}$, $\partial_{LE}$, $\partial_{T_s}$, $\partial_{T_c}$, and $\partial_{T_{ao}}$ were more sensitive to the $\theta$ change, but less sensitive to the LAI change when LAI > 1.0.

Fig. S3 shows the global patterns of $\partial_{\emptyset_n}$, $\partial_{LE}$, $\partial_{T_s}$, $\partial_{T_c}$, and $\partial_{T_{ao}}$. The $T_{af}$ was more positively sensitive to $\emptyset_n$ in the November–January seasonal period (0.011 ± 0.01) than in the May–July seasonal period (0.010 ± 0.008) (Fig. S3 a1 and a2), while $T_{af}$ was less positively sensitive to $T_{ao}$ in the November–January period (2.51 ± 0.05) than in the May–July period (3.03 ± 0.05) (Fig. S3 e1 and e2). LE, $T_s$, and $T_c$ showed less negative sensitivity to $T_{af}$ in the November–January period ($\partial_{LE}$: -0.011 ± 0.0008, $\partial_{T_s}$: -0.087 ± 0.0008, $\partial_{T_c}$: -1.495 ± 0.06) than in the May–July period ($\partial_{LE}$: -0.012 ± 0.0006, $\partial_{T_s}$: -0.088 ± 0.0009, $\partial_{T_c}$: -1.961 ± 0.09) (Fig. S3 b1, b2, c1, c2, d1, and d2). In general, the absolute values of $\partial_{\emptyset_n}$, $\partial_{LE}$, $\partial_{T_s}$, $\partial_{T_c}$, and $\partial_{T_{ao}}$ were higher at the northern and southern high latitudes than those of the low tropical latitudes (Fig. S3), indicating a greater sensitivity of high-latitude $T_{af}$ to the environmental changes due to a greater $\theta$ (Fig. 6). $\partial_{\emptyset_n}$, $\partial_{LE}$, $\partial_{T_s}$, $\partial_{T_c}$, and $\partial_{T_{ao}}$ at higher latitudes showed greater interannual variability than at low latitudes (Fig. S3). This could be one of the main reasons why average $\Delta T_a$ at northern/southern high latitudes has greater seasonal dynamics than $\Delta T_a$ at low tropical latitudes (Fig. 7a).

The above principle analysis suggests that the seasonality of two non-climatic factors – solar geometry (i.e., $\theta$) and canopy phenology (i.e., LAI) – probably controls the seasonal variabilities of $\Delta T_a$. As shown in Fig. 7, both $\theta$ and LAI exhibited highly consistent seasonality with $\Delta T_a$ at northern and southern high latitudes. In the May–July period, the forests at northern high latitudes experienced the lowest $\theta$ and highest canopy LAI (Fig. 7b and c). These tended to induce the most notable differences in energy fluxes between forest ecosystems and open air spaces, resulting in the highest $\Delta T_a$ (Fig. 7a). Given the inverse seasonality in the southern hemisphere, the forests at northern high latitudes had the highest $\theta$ and lowest LAI (Fig. 7b and c), exhibiting the lowest $\Delta T_a$ between forests and open air spaces (Fig. 7a); the opposite occurred during the November–January period for both the northern and southern high latitudes. The low-latitude tropical forests revealed little seasonality of $\theta$ and LAI, as well as of $\Delta T_a$ (Fig. 7).

## 3.4. $T_c$ and $T_{ao}$ are the two most important climatic factors in impacting $T_{af}$ and $\Delta T_a$

Fig. 8 shows the global patterns of $C_{\emptyset_n}$, $C_{LE}$, $C_{T_s}$, $C_{T_c}$, and $C_{T_{ao}}$, which follow the order of: $C_{T_{ao}}$ (51.8–68.3%) > $C_{T_c}$ (30.70–42.60%) > $C_{T_s}$ (0.27–8.70%) > $C_{\emptyset_n}$ (0.22–2.08%) > $C_{LE}$ (0.00–0.13%) (Fig. S4). $C_{T_{ao}}$ and $C_{T_s}$ are the two most important climate factors in mediating $T_{af}$. At the northern high latitudes, the minimum $C_{T_{ao}}$ (58.23 ± 1.88%) and

<caption>**Table 2**</caption>

Proportions of simulated $T_{af}$ located at different thresholds in the May–July and November–January periods at the northern high latitudes (> 23.5°N), tropics (23.5°S–23.5°N) and southern high latitudes (> 23.5°S).

<table><thead><tr><th rowspan="2">Range of T<sub>af</sub> (°C)</th><th colspan="2">Northern high latitudes</th><th colspan="2">Tropics</th><th colspan="2">Southern high latitudes</th></tr><tr><th>May–Jul</th><th>Nov–Jan</th><th>May–Jul</th><th>Nov–Jan</th><th>May–Jul</th><th>Nov–Jan</th></tr></thead><tbody><tr><td>&lt; 0.0</td><td>8.0%</td><td>66.5%</td><td>0.2%</td><td>1.1%</td><td>2.3%</td><td>1.2%</td></tr><tr><td>0–10</td><td>58.5%</td><td>27.0%</td><td>5.2%</td><td>3.3%</td><td>74.0%</td><td>10.3%</td></tr><tr><td>10–20</td><td>29.8%</td><td>6.4%</td><td>57.3%</td><td>13.3%</td><td>23.7%</td><td>41.0%</td></tr><tr><td>&gt; 20</td><td>3.7%</td><td>0.2%</td><td>37.3%</td><td>82.3%</td><td>0.0%</td><td>47.5%</td></tr></tbody></table>
<caption>**Table 3**</caption>

Statistics of simulated ΔT<sub>a</sub> located at different thresholds in the May–July and November–January periods at the northern high latitudes (> 23.5°N), tropics (23.5°S–23.5°N) and southern high latitudes (> 23.5°S).

<table><thead><tr><th rowspan="2">Range of ΔT<sub>a</sub> (°C)</th><th colspan="2">Northern high latitudes</th><th colspan="2">Tropics</th><th colspan="2">Southern high latitudes</th></tr><tr><th>May–Jul</th><th>Nov–Jan</th><th>May–Jul</th><th>Nov–Jan</th><th>May–Jul</th><th>Nov–Jan</th></tr></thead><tbody><tr><td>&lt; -5.0</td><td>51.8%</td><td>8.5%</td><td>20.2%</td><td>21.6%</td><td>1.7%</td><td>29.8%</td></tr><tr><td>-5 to -2.5</td><td>43.1%</td><td>7.8%</td><td>38.9%</td><td>19.5%</td><td>32.7%</td><td>20.2%</td></tr><tr><td>-2.5–0.0</td><td>4.4%</td><td>22.0%</td><td>38.4%</td><td>20.6%</td><td>50.7%</td><td>30.3%</td></tr><tr><td>&gt; 0.0</td><td>0.7%</td><td>61.7%</td><td>2.5%</td><td>38.3%</td><td>14.9%</td><td>19.7%</td></tr><tr><td>Mean</td><td>-5.95</td><td>-0.24</td><td>-3.52</td><td>-3.83</td><td>-1.64</td><td>-6.84</td></tr><tr><td>Median</td><td>-5.20</td><td>0.66</td><td>-2.87</td><td>-1.74</td><td>-1.94</td><td>-5.93</td></tr><tr><td>Standard deviation (Std)</td><td>3.15</td><td>3.74</td><td>2.49</td><td>6.17</td><td>1.78</td><td>6.62</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Sensitivity of $T_{af}$ to $\emptyset_n(\partial_{\emptyset_n})$, $LE(\partial_{LE})$, $T_s(\partial_{T_s})$, $T_c(\partial_{T_c})$ and $T_{ao}(\partial_{T_{ao}})$ under different LAIs.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 6. Contributions of $\bigcirc_n(C_{\bigcirc_n})$, $LE(C_{LE})$, $T_s(C_{T_s})$, $T_c(C_{T_c})$ and $T_{ao}(C_{T_{ao}})$ to $T_{af}$ under different $\theta$.</figcaption></figure>

maximum $C_{T_c}$ (40.49 ± 1.83%) occurred in the May–July period, while the maximum $C_{T_{ao}}$ (62.14 ± 1.64%) and minimum $C_{T_s}$ (35.97 ± 1.37%) occurred in the November–January period. The $C_{T_{ao}}$ and $C_{T_c}$ at southern high latitudes showed inverse interannual variability compared to those at the northern high latitudes. The $C_{T_{ao}}$ and $C_{T_c}$ of tropical forests at low latitudes showed little seasonal variability. In addition, tropical forests, in general, had higher $C_{T_{ao}}$ and lower $C_{T_c}$ than those at high latitudes throughout the whole year.

The aforementioned analyses revealed that the $T_c$ and $T_{ao}$ terms, which are deconstructed from sensible heat fluxes ($H_{soil→af}$, $H_{af→canopy}$ and $H_{canopy→ao}$) and canopy longwave radiative heat fluxes

($LW_{canopy}$ [$1 - \exp(-\frac{cLAI}{u})$]), are two dominant contributors in controlling $T_{af}$. This suggests that $T_c$ and $T_{ao}$ are likely to be important climatic factors that mediate the direction and magnitude of $\Delta T_a$. To demonstrate this, we plotted the differences between $T_c$ and $T_{ao}$ (i.e., $T_c - T_{ao}$) against corresponding $\Delta T_a$. Fig. S5 is the scatter plot between $T_c - T_{ao}$ and field $\Delta T_a$ from 373 published observations; Fig. 9, instead, shows the global simulated results. It is worth noting that $T_c - T_{ao}$ was well correlated with the corresponding $\Delta T_a$ from both field observations ($R^2 = 0.51, p < 0.001, N = 373$) and global simulated results (May–July: $R^2 = 0.35, p < 0.001, N = 2843$; November–January: $R^2 = 0.51$,
<figure><img src="image_9.png" /><figcaption>Fig. 7. Seasonality of ΔT_a, LAI and θ at the northern high latitudes (> 23.5°N), Tropics (23.5°S–23.5°N) and southern high latitudes (> 23.5°S). (a) ΔT_a; (b) LAI; (c) θ.</figcaption></figure>

p < 0.001, N = 2223). Forests in the May–July period revealed higher $T_c - T_{ao}$ (> 99.5% positive forest lands) and resulted in positive ΔT_a when $T_c - T_{ao} < -10.0$ °C; however, in the November–January period 87.7% of the forests exhibited negative $T_c - T_{ao}$ and resulted in positive ΔT_a when $T_c - T_{ao} < -0.0$ °C.

## 4. Discussion

There has been growing interest in the impact of forest ecosystem on human health (Wolff et al., 2018), as forests not only impact the Earth's climate at large scale but also play an important role in regulating the local temperatures (Ellison et al., 2017). The loss of forest cooling services, via deforestation, poses a health risk to tropical rural villagers through increased heat exposure (Wolff et al., 2018; Masuda et al., 2019; Suter et al., 2019). The results of our study confirm the claims of the published literature we reviewed. Our findings show that regional/ global investigations based on satellite-observed land surface

temperatures (i.e., MODIS), which is measured at tree canopy, likely overestimate the air temperature experienced by humans by considerable amount. The impact of this overestimation can be significant, as forest canopies cover > 30% of the global land surface (FAO, 2016), not to mention the hundreds of thousands of forest spaces in urban areas (Chen et al., 2012). For example, field observations in Shenzhen's urban park in subtropical southern China (longitude: 113.97°E, latitude: 22.57°N) revealed that in hot summer weather ($T_{af}$ = 32.25 °C) forests are 5.60 °C cooler than open air spaces approximately 20% of the time (Lei et al., 2011). Clearly, the biophysical effects of forests can, to some extent, reduce the magnitude of $T_{af}$ fluctuation caused by background climate changes. Preserving tree cover, especially in tropical/subtropical landscapes, can reduce high air temperatures and buffer some of the extremes that are otherwise likely to arise with climate change (Ellison et al., 2017). Previous satellite-based observations, land-atmosphere model simulations and field investigations have attributed the inconsistent forest cooling/warming at temperate and boreal latitudes to the complex role of albedo (Bonan, 2008; Zeng et al., 2017). However, many of these studies overlooked the diverse temperature variables that were considered; for example, satellite-retrieved land surface temperatures or $T_s$ products (Li et al., 2015), modeled near-surface temperature (Lee et al., 2011), and field observed air temperature (Chang et al., 2007). It is necessary for studies to adjust satellite-derived $T_s$ and model-derived near-surface temperatures to more accurately reflect real air temperature that human experiences.

Note also that our study reveals that the biophysical effects of forests on air temperature and soil surface temperature can be quite disparate. Numerous studies have reported that ET (Li et al., 2015), aerodynamic roughness (Lee et al., 2011; Burakowski et al., 2018), albedo (Lee et al., 2011), surface resistance (Li et al., 2019), or background climate (Pitman et al., 2011) are several major factors that affect the forest cooling/warming on large-scale, near-surface temperatures and land surface temperatures (Taha et al., 1991; Schwartz, 1996; Pearlmutter et al., 1999; Fitzjarrald et al., 2001; Shashua-Bar and Hoffman, 2002; Dimoudi and Nikolopoulou, 2003; Jonsson, 2004; Bowler et al., 2010; Zhao et al., 2014). Our study suggests that the $T_c$ and $T_{ao}$ terms are probably the two dominant climatic factors influencing the biophysical effects of forests on local air temperature. The solar geometry (i.e., θ) and canopy phenology (i.e., LAI) can also influence the energy that penetrates the forest canopy, redistributing the energy between the canopy and understory layers (Holmes and Smith, 1977; Jackson et al., 1981; Running et al., 1989; Maia and Loureiro, 2005; Li et al., 2015) and, subsequently, mediating the biophysical effects of forests on seasonal ΔT_a. This is an important outcome of our analysis, which depicts the diverse impacts of the biophysical effects of forests on local air temperature rather than soil surface temperature.

## 5. Conclusions

By partitioning forest lands into three vertical layers, the CAS land surface energy balance model successfully decomposes the sensible heat fluxes and longwave radiations within the tree canopy, forest air space, and soil surface layers of forest ecosystems. The improvement of energy balance makes it possible to separate the individual biophysical effects of forests on local air temperatures. Global applications of the novel CAS model suggest that forests at northern and southern high latitudes exhibit larger quantities of negative ΔT_a (i.e., cooling effects) in warm seasons and positive ΔT_a (i.e., warming effects) in cool seasons compared with the low-latitude forests in tropical/subtropical regions showing less seasonal dynamics. Furthermore, contribution analysis suggests that the differences between $T_c$ and $T_{ao}$, two dominant variables of climate to $T_{af}$, explain a large part of forest cooling/warming (May–July: R² = 0.35; November–January: R² = 0.51), and their seasonal variabilities are jointly regulated by both the seasonality of solar geometry (i.e., θ) and canopy phenology (i.e., LAI).

In summary, this study proposes a framework to decompose the
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Contributions of $\emptyset_n(C_{\emptyset_n})$, $LE(C_{LE})$, $T_s(C_T)$, $T_c(C_T)$ and $T_{ao}(C_{T_{ao}})$ to $T_{af}$ at the northern high latitudes (> 23.5°N), tropics (23.5°S–23.5°N) and southern high latitudes (> 23.5°S).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Scatter plots between $T_c - T_{ao}$ and $\Delta T_{ao}$. (a) May–July period; (b) November–January period.</figcaption></figure>

biophysical effects of forests on local air temperatures using a novel three-layered CAS model with robust performance. It is anticipated that the CAS model can be widely adopted for global applications. The collected site observations are mostly located in urban areas and are skewed towards northern high latitudes; and more in-situ data from natural forests are needed for intensive model calibrations in the future.

**Declaration of competing interest**

The authors declare no conflict of interest.

**Acknowledgements**

We thank the editor and the anonymous reviewers for their insights and valuable comments that helped improve the quality of the manuscript. This study was supported by the National Natural Science Foundation of China [grant numbers 41971275, 31971458, 41430529], Special High-level Plan Project of Guangdong Province [grant number 2016TQ03Z354], the Key Programs of the Chinese Academy of Sciences [grant numbers KFZD-SW-312, QYZDJ-SSW-DQC003], the Pearl River S &T Nova Program of Guangzhou [grant number 201610010134], the GDAS' project of Science and Technology Development [grant numbers 2016GDASRC-0101, 2017GDASCX-0501, 2019GDASYL-0103002], and
the Water Resource Science and Technology Innovation Program of
Guangdong Province [grant number 2016-16].

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://
doi.org/10.1016/j.envint.2019.105080.

References

Alkama, R., Cescatti, A., 2016. Biophysical climate impacts of recent changes in global forest cover. Science 351 (6273), 600–604.

Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998. Crop evapotranspiration-guidelines for computing crop water requirements-FAO irrigation and drainage paper 56. FAO, Rome 300 (9), D05109.

Anderson, R.G., Canadell, J.G., Randerson, J.T., Jackson, R.B., Hungate, B.A., Baldocchi, D.D., Diffenbaugh, N.S., 2011. Biophysical considerations in forestry for climate protection. Front. Ecol. Environ. 9 (3), 174–182.

Berrisford, P., Källberg, P., Kobayashi, S., Dee, D., Uppala, S., Simmons, A.J., ... Sato, H., 2011. Atmospheric conservation properties in ERA-Interim. Quarterly Journal of the Royal Meteorological Society 137 (659), 1381–1399.

Bonan, G.B., 2008. Forests and climate change: forcings, feedbacks, and the climate benefits of forests. Science 320 (5882), 1444–1449.

Bounoua, L., Defries, R., Collatz, G.J., Sellers, P., Khan, H., 2002. Effects of land cover conversion on surface climate. Clim. Chang. 52 (1–2), 29–64.

Bowler, D.E., Buyung-Ali, L., Knight, T.M., Pullin, A.S., 2010. Urban greening to cool towns and cities: a systematic review of the empirical evidence. Landsc. Urban Plan. 97, 147–155.

Burakowski, E., Tawfik, A., Ouimette, A., Lepine, L., Novick, K., Ollinger, S., Bonan, G., 2018. The role of surface roughness, albedo, and Bowen ratio on ecosystem energy balance in the Eastern United States. Agric. For. Meteorol. 249, 367–376.

Ca, V.T., Asaeda, T., Abu, E.M., 1998. Reductions in air conditioning energy caused by a nearby park. Energy Build. 29 (1), 83–92.

Cao, X., Onishi, A., Chen, J., Imura, H., 2010. Quantifying the cool island intensity of urban parks using ASTER and IKONOS data. Landsc. Urban Plan. 96 (4), 224–231.

Chang, C.R., Li, M.H., Chang, S.D., 2007. A preliminary study on the local cool-island intensity of Taipei city parks. Landsc. Urban Plan. 80, 386–395.

Chen, X.Z., Su, Y.X., Li, D., Huang, G.Q., Chen, W.Q., Chen, S.S., 2012. Study on the cooling effects of urban parks on surrounding environments using Landsat TM data: a case study in Guangzhou, southern China. Int. J. Remote Sens. 33 (18), 5889–5914.

Choudhury, B.J., Monteith, J.L., 1988. A four-layer model for the heat budget of homo-
geneous land surfaces. Quarterly Journal of the Royal Meteorological Society 114
(480), 373–398.

Devaraju, N., de Noblet-Ducoudré, N., Quesada, B., Bala, G., 2018. Quantifying the relative importance of direct and indirect biophysical effects of deforestation on surface temperature and teleconnections. J. Clim. 31 (10), 3811–3829.

Dimoudi, A., Nikolopoulou, M., 2003. Vegetation in the urban environment: microclimatic analysis and benefits. Energy Build. 35, 69–76.

Dooge, J.C., 1992. Sensitivity of runoff to climate change: a Hortonian approach. Bull. Am. Meteorol. Soc. 73 (12), 2013–2024.

Ellison, D., Morris, C.E., Locatelli, B., Sheil, D., Cohen, J., Murdiyarso, D., Gaveau, D., 2017. Trees, forests and water: cool insights for a hot world. Glob. Environ. Chang. 43, 51–61.

Feyisa, G.L., Dons, K., Milby, H., 2014. Efficiency of parks in mitigating urban heat island effect: an example from Addis Ababa. Landsc. Urban Plan. 123, 87–95.

Fitzjarrald, D.R., Acevedo, O.C., Moore, K.E., 2001. Climatic consequences of leaf presence in the eastern United States. J. Clim. 14 (4), 598–614.

Food and Agricultural Organization (FAO), 2016. Global Forest Resources Assessment 2015: How Are the World's Forests Changing? FAO, Rome.

Friedlingstein, P., Dufresne, J.L., Cox, P.M., Rayner, P., 2003. How positive is the feedback between climate change and the carbon cycle? Tellus B Chem. Phys. Meteorol. 55 (2), 692–700.

Georgi, N.J., Zafiradias, K., 2006. The impact of park trees on microclimate in urban areas. Urban Ecosyst. 9, 195–209.

He, L., Chen, J.M., Liu, J., Mo, G., Joiner, J., 2017. Angular normalization of gome-2 sun-induced chlorophyll fluorescence observation as a better proxy of vegetation productivity. Geophys. Res. Lett. 5691–5699.

Heidinger, A.K., Foster, M.J., Walther, A., Zhao, X., 2014. The pathfinder atmospheres-extended AVHRR climate dataset. Bulletin of the American Meteorological Society 95 (6), 909–922.

Holmes, M.G., Smith, H., 1977. The function of phytochrome in the natural environment—II. The influence of vegetation canopies on the spectral energy distribution of natural daylight. Photochem. Photobiol. 25 (6), 539–545.

Jackson, R.D., Idso, S.B., Reginato, R.J., Pinter, P.J., 1981. Canopy temperature as a crop water stress indicator. Water Resour. Res. 17 (4), 1133–1138.

Jauregui, E., 1991. Influence of a large urban park on temperature and convective precipitation in a tropical island. Energy Build. 15 (3), 457–463.

Ji, P., Zhu, C.Y., Li, S.H., 2012. Effects of urban river width on the temperature and humidity of nearby green belts in summer. Chin. J. Appl. Ecol. 23 (3), 679.

Jones, H.G., Rotenberg, E., 2011. Energy, Radiation and Temperature Regulation in Plants. eLS.

Jonsson, P., 2004. Vegetation as an urban climate control in the subtropical city of Gaborone, Botswana. Int. J. Climatol. 24, 1307–1322.

Lee, X., Goulden, M.L., Hollinger, D.Y., Barr, A., Black, T.A., Bohrer, G., ... Katul, G., 2011. Observed increase in local cooling effect of deforestation at higher latitudes. Nature 479 (7373), 384–387.

Lei, J.L., Liu, T., Yan-Yan, W.U., Zhuang, X.Y., Xie, L.S., 2011. Effects of structure characteristics of urban green land on the temperature-lowering in Shenzhen City. J. Northwest For. Univ. 26 (4), 218–223.

Li, Y., Zhao, M., Motesharrei, S., Mu, Q., Kalnay, E., Li, S., 2015. Local cooling and warming effects of forests based on satellite observations. Nat. Commun. 6, 6603.

Li, D., Liao, W., Rigden, A.J., Liu, X., Wang, D., Malyshev, S., Shevliakova, E., 2019. Urban heat island: aerodynamics or imperviousness? Sci. Adv. 5 (4), eaau4299.

Liu, J.M., Li, S.H., Yang, Z.F., 2008. Temperature and humidity effect of urban green spaces in Beijing in summer. Chin. J. Ecol. 27 (11), 1972–1978.

Liu, L., Liao, J., Chen, X., Zhou, G., Su, Y., Xiang, Z., Xiong, X., 2017. The Microwave Temperature Vegetation Drought Index (MTVDI) based on AMSR-E brightness temperatures for long-term drought assessment across China (2003–2010). Remote Sens. Environ. 199, 302–320.

Maia, A.S.C., Loureiro, C.B., 2005. Sensible and latent heat loss from the body surface of Holstein cows in a tropical environment. Int. J. Biometeorol. 50 (1), 17–22.

Marland, G., Pielke Sr., R.A., Apps, M., Avissar, R., Betts, R.A., Davis, K.J., Katzenberger, J., 2003. The climatic impacts of land surface change and carbon management, and the implications for climate-change mitigation policy. Clim. Pol. 3 (2), 149–157.

Masuda, Y.J., Castro, B., Aggraeni, I., Wolff, N.H., Ebi, K., Garg, T., Game, E.T., Krenz, J., Spector, J., 2019. How are healthy, working populations affected by increasing temperatures in the tropics? Implications for climate change adaptation policies. Glob. Environ. Chang. 56, 29–40.

Monteith, J., Unsworth, M., 2007. Principles of Environmental Physics. Academic Press.

Mu, Q., Zhao, M., Running, S.W., 2011. Improvements to a MODIS global terrestrial evapotranspiration algorithm. Remote Sensing of Environment 115 (8), 1781–1800.

Myneni, R.B., Hoffman, S., Knyazikhin, Y., Privette, J.L., Glassy, J., Tian, Y., et al., 2002. Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data. Remote Sensing of Environment 83 (1–2), 214–231.

Norman, J.M., Kustas, W.P., Humes, K.S., 1995. Source approach for estimating soil and vegetation energy fluxes in observations of directional radiometric surface temperature. Agric. For. Meteorol. 77 (3–4), 263–293.

Oke, T.R., 1989. The micrometeorology of the urban forest. Philos. Trans. R. Soc. B 324, 335–349.

Oliveira, S., Andrade, H., Vaz, T., 2011. The cooling effect of green spaces as a contribution to the mitigation of urban heat: a case study in Lisbon. Build. Environ. 46, 2186–2194.

Paltridge, G.W., 1980. Cloud-radiation feedback to climate. Q. J. R. Meteorol. Soc. 106 (450), 895–899.

Pearlmutter, D., Bitan, A., Berliner, P., 1999. Microclimatic analysis of “compact”urban canyons in an arid zone. Atmos. Environ. 33, 4143–4150.

Peng, S.S., Piao, S., Zeng, Z., Ciais, P., Zhou, L., Li, L.Z., Myneni, R.B., Yin, Y., Zeng, H., 2014. Afforestation in China cools local land surface temperature. Proc. Natl. Acad. Sci. 111 (8), 2915–2919.

Pitman, A.J., Avila, F.B., Abramowitz, G., Wang, Y.P., Phipps, S.J., de Noblet-Ducoudré, N., 2011. Importance of background climate in determining impact of land-cover change on regional climate. Nat. Clim. Chang. 1 (9), 472.

Potchter, O., Cohen, P., Bitan, A., 2006. Climatic behavior of various urban parks during a hot and humid summer in the Mediterranean city of Tel Aviv, Israel. Int. J. Climatol. 26, 1695–1711.

Potchter, O., Goldman, D., Iluz, D., Kadish, D., 2012. The climatic effect of a manmade oasis during winter season in a hyper arid zone: the case of southern Israel. J. Arid Environ. 87 (12), 231–242.

Running, S.W., Nemani, R.R., Peterson, D.L., Band, L.E., Potts, D.F., Pierce, L.L., Spanner, M.A., 1989. Mapping regional forest evapotranspiration and photosynthesis by coupling satellite data with ecosystem simulation. Ecology 70 (4), 1090–1101.

Schwartz, M.D., 1996. Examining the spring discontinuity in daily temperature ranges. J. Clim. 9 (4), 803–808.

Shashua-Bar, L., Hoffman, M.E., 2002. The green CTTC model for predicting the air temperature in small urban wooded sites. Build. Environ. 37, 1279–1288.

Shashua-Bar, L., Pearlmutter, D., Erell, E., 2009. The cooling efficiency of urban landscape strategies in a hot dry climate. Landsc. Urban Plan. 92, 179–186.

Shashua-Bar, L., Tsiros, I.X., Hoffman, M.E., 2010. A modeling study for evaluating passive cooling scenarios in urban streets with trees. Case study: Athens, Greece. Build. Environ. 45, 2798–2807.

Simard, M., Pinto, N., Fisher, J.B., Baccini, A., 2011. Mapping forest canopy height globally with spaceborne lidar. J. Geophys. Res. 116https://doi.org/10.1029/2011JG001708. G04021.

Susca, T., Gaffin, S.R., Dell'Osso, G.R., 2011. Positive effects of vegetation: urban heat island and green roofs. Environ. Pollut. 159, 2119–2126.

Suter, M., Miller, K.A., Anggraeni, I., Ebi, K.L., Game, E.E., Krenz, J., Masuda, Y.J., Sheppard, L., Wolff, N.H., Spector, J., 2019. Association between work in deforested, compared to forested, areas and human heat strain: an experimental study in a rural tropical environment. Environ. Res. Lett. 14, 084012.

Taha, H., Akbari, H., Rosenfeld, A., 1991. Heat island and oasis effects of vegetative canopies: micro-meteorological field-measurements. Theor. Appl. Climatol. 44 (4), 123–138.

Vidrih, B., Medved, S., 2013. Multiparametric model of urban park cooling island. Urban For. Urban Green. 12 (2), 220–229.

Wei, Y., Liu, S., Huntzinger, D.N., et al., 2014. The North American Carbon Program Multi-scale Synthesis and Terrestrial Model Intercomparison Project - Part 2: Environmental driver data[J]. Geoscientific Model Development Discussions 6 (4), 5375–5422.

Wan, Z., 2014. New refinements and validation of the modis land-surface temperature/
emissivity products. Remote Sensing of Environment 140 (1), 36–45.

Wolff, N.H., Masuda, Y.J., Meijaard, E., Wells, J.A., Game, E.T., 2018. Impacts of tropical deforestation on local temperature and human well-being perceptions. Glob. Environ. Chang. 52, 181–189.

Wong, N.H., Yu, C., 2005. Study of green areas and urban heat island in a tropical city. Habitat Int. 29 (3), 547–558.

Wu, C.L., Chau, K.W., Huang, J.S., 2007. Modelling coupled water and heat transport in a soil–mulch–plant–atmosphere continuum (SMPAC) system. Appl. Math. Model. 31 (2), 152–169.

Zeng, Z., Piao, S., Li, L.Z., Zhou, L., Ciais, P., Wang, T., Mao, J., 2017. Climate mitigation from vegetation biophysical feedbacks during the past three decades. Nat. Clim. Chang. 7 (6), 432.

Zhang, Z., Lv, Y., Pan, H., 2013. Cooling and humidifying effect of plant communities in subtropical urban parks. Urban For. Urban Green. 12 (3), 323–329.

Zhang, M., Lee, X., Yu, G., Han, S., Wang, H., Yan, J., Kim, J., 2014. Response of surface air temperature to small-scale land clearing across latitudes. Environ. Res. Lett. 9 (3), 034002.

Zhao, L., Lee, X., Smith, R.B., Oleson, K., 2014. Strong contributions of local background climate to urban heat islands. Nature 511 (7508), 216–219.

Zhou, G., Wei, X., Chen, X., Zhou, P., Liu, X., Xiao, Y., Su, Y., 2015. Global pattern for the effect of climate and land cover on water yield. Nat. Commun. 6, 5918.

Zhou, P., Li, Q., Zhou, G., Wei, X., Zhang, M., Liu, Z., Scott, D.F., 2018. Reply to ‘Flawed assumptions compromise water yield assessment’. Nat. Commun. 9 (1), 4788.