Dynamic Cooling Effects of Permanent Urban Green Spaces in Beijing, China

Urban green spaces (UGSs) play a critical role in human thermal comfort, energy consumption and urban ecology. Although the heat mitigation capability of UGSs has been frequently reported, many of the current understandings are based on short-term observations, and the long-term temporal dynamics of UGS cooling effects are still lacking. This gap may cause over- or underestimation and largely ignores how the cooling effects change with climate change and urban growth. Accordingly, we used Landsat-based time series data to analyze the changes in permanent UGS greenness, surface-cooling effects and their biophysical responses in Beijing in the past 40 years (1984–2020). The results demonstrate segmented changes in UGS surface cooling that were mainly linked to the responses of canopy transpiration and albedo to vegetation conditions. During a rapid greening of UGSs in the recent two decades, transpiration cooling dominated albedo-induced warming to provide a discernable cooling enhancement. In addition, such enhancement showed seasonal differences ranging from less than 1 °C to more than 2 °C, and the most evident enhancement occurred on summer days (~2.4 °C) when vegetation is most needed to provide cooling. The highlighted dynamics of UGSs help urban planners better balance the maintenance costs and the environmental gains for UGS management.


Introduction
Urban green spaces (UGSs) are defined as vegetated open spaces in urban areas [1]. They improve environmental sustainability and residents' well-being in multiple ways. For instance, UGSs play an essential role in maintaining high levels of urban biodiversity [2]. They also reduce air pollutants and mitigate urban floods [3,4]. Additionally, UGSs facilitate people's mental health and reduce mortality by providing elevated exposure to nature, which helps alleviate stress and enhance physical activity [5,6]. Another major function of UGSs is climate mitigation. It is well known that UGSs reduce urban heat through both shading and transpiration. A large tree canopy with dense leaves forms solid shadows that prevent heating from solar radiation, while the process of transpiration consumes this solar energy to evaporate liquid water from leaves to the atmosphere, consequently reducing both the leaf surface temperature and surrounding air temperature [7]. This cooling effect of UGSs will be increasingly important in future urban development, since urban areas have been projected to experience more serious warming than regional background warming [8].
There are currently three main methods used to measure UGS cooling effects. First, tree-level field measurements are conducted for various species to find potential drivers (such as leaf area index (LAI), canopy geometry, stress tolerance, microclimate factors and growing conditions) that affect the cooling efficiency [9][10][11]. Second, remotely sensed images at the city scale are widely utilized to map UGSs [12] and explore the relationship between UGS cooling intensity/distance and their bio and geometric features [13][14][15]. The cooling intensity is defined as the thermal difference between the UGS and its surroundings [16], while the cooling distance measures the spatial extent beyond which the UGS ceases to provide local cooling effects [17]. Third, the cooling effects at the neighborhood scale are estimated by mechanistic models, which consider street scenarios with different planting patterns and urban morphologies [18,19]. Despite the remarkable contributions of these studies, most of them were conducted during a single period (from a few months [20] to a few years [11]), and the long-term temporal dynamics of permanent UGS cooling effects have rarely been reported. In drylands, plant physiology is tightly linked to changes in precipitation and frequency [21,22]. For plants in dry cities, their physiological traits could be more variable due to the need for human management and irrigation [23]. Consequently, changes in vegetation greenness (e.g., greening or browning) alter land-atmosphere interactions through biophysical processes, and the cooling effects are therefore modified over time. For instance, a field campaign conducted in a semiarid shrub ecosystem observed that deep-soil moisture led to an increase in LAI and a decrease in canopy albedo [24]. In this process, a higher LAI increases canopy transpiration to result in more cooling effects, while a lower albedo absorbs more solar energy to cause warming effects. Finally, the two biophysical responses collectively determine the change in vegetation thermal effects [25].
Remote sensing techniques intrinsically have advantages in land surface monitoring due to satellite repeatability. Landsat sensors have circled the Earth since 1972, providing the longest and richest continuous satellite observations. With a spatial resolution of 30 m, Landsat images are capable of featuring finer-scale anthropogenic impacts and detailed UGS (e.g., parks, urban forests, street trees, lawns) characteristics and are thus widely used for UGS monitoring and modeling [15,26,27]. Landsat time series data provide more dynamic and comprehensive land processes than bitemporal image comparisons and capture abrupt, stochastic events, as well as subtle changes, over long periods [28]. Hence, the data are well suited for the long-term monitoring of UGSs and can be used to report how UGSs grow, decline and potentially change with urbanization and management.
Based on Landsat-based observations, this study performed time series analyses (1984-2020) on permanent UGSs in Beijing, a typical semiarid metropolis located in northern China. Seasonal cycles and stochastic noise were removed to extract the trends of various surface parameters, including the vegetation index, land surface temperature (LST), evapotranspiration (ET) and albedo. We also applied segmented regressions to time series data to find potential breakpoints. In general, we aimed to explore the dynamics of permanent UGSs and the changes in their environmental benefits by answering the following questions: (i) did the greenness of these UGSs change over the past 40 years, and when did the changes happen?; (ii) did these changes provide variant cooling effects?; and (iii) how do vegetation physiological processes account for these changing cooling effects?

Study Area
Beijing is one of the leading metropolises in China and has over 3000 years of history. As of 2018, there were more than 21 million people who resided in Beijing. The city is located on the North China Plain between 39 • 28 -41 • 05 N and 115 • 25 -117 • 30 E, experiencing a warm temperate semiarid continental monsoon climate that features four distinct seasons. The city has an annually averaged daily temperature of 10-12 • C, with monthly mean values varying between −5 • C and 30 • C. The annual precipitation in Beijing is approximately 650 mm, which mostly occurs in July and August. The growing season typically starts in April and lasts until the end of October.
The study was carried out in permanent UGSs in Beijing located within the 4th Ring Road (Figure 1). These UGSs were identified from the Chinese urban park dataset by Li et al. [29] and the Google historical images. We selected the UGSs that cover at least 100 Landsat pixels (approximately equal to 10 ha) because too small UGSs are easily affected by their surroundings and thus raise uncertainty. Formed before 1949, most of these UGSs are famous historical parks, e.g., North Lake Park and Coal Hill Park, consisting of abundant arbors, shrubs, herbaceous species and water bodies [30]. Due to their outstanding historical and cultural values, both landscape patterns and garden elements are well protected. Thus, these permanent UGSs provide us with a suitable testbed to quantify the long-term cooling dynamics of urban vegetation.
Beijing is one of the leading metropolises in China and has over 3000 years of histor As of 2018, there were more than 21 million people who resided in Beijing. The city located on the North China Plain between 39°28′-41°05′ N and 115°25′-117°30′ E, exper encing a warm temperate semiarid continental monsoon climate that features four distin seasons. The city has an annually averaged daily temperature of 10-12 °C, with month mean values varying between −5 °C and 30 °C. The annual precipitation in Beijing is a proximately 650 mm, which mostly occurs in July and August. The growing season typ cally starts in April and lasts until the end of October.
The study was carried out in permanent UGSs in Beijing located within the 4th Rin Road (Figure 1). These UGSs were identified from the Chinese urban park dataset by et al. [29] and the Google historical images. We selected the UGSs that cover at least 10 Landsat pixels (approximately equal to 10 ha) because too small UGSs are easily affecte by their surroundings and thus raise uncertainty. Formed before 1949, most of these UGS are famous historical parks, e.g., North Lake Park and Coal Hill Park, consisting of abu dant arbors, shrubs, herbaceous species and water bodies [30]. Due to their outstandin historical and cultural values, both landscape patterns and garden elements are well pr tected. Thus, these permanent UGSs provide us with a suitable testbed to quantify th long-term cooling dynamics of urban vegetation.  Figure 2 illustrates the overall workflow of the study. We first combined Landsat and Landsat 8 surface reflectance datasets to form a long-term time series from 1984  Figure 2 illustrates the overall workflow of the study. We first combined Landsat 5 and Landsat 8 surface reflectance datasets to form a long-term time series from 1984 to 2020. Next, based on the dataset, we retrieved four surface parameters including the normalized difference vegetation index (NDVI), LST, ET and albedo using different algorithms. Finally, we analyzed the greenness and cooling dynamics of UGSs and their biophysical responses using time series tools. Detailed descriptions of these steps are shown in the following sections. 2020. Next, based on the dataset, we retrieved four surface parameters including the normalized difference vegetation index (NDVI), LST, ET and albedo using different algorithms. Finally, we analyzed the greenness and cooling dynamics of UGSs and their biophysical responses using time series tools. Detailed descriptions of these steps are shown in the following sections.

Figure 2.
Overall workflow of the study. The main steps consist of data preprocessing, parameter retrieval and dynamic cooling assessment of UGSs. Note that SR refers to surface reflectance.

Landsat Surface Reflectance Collection
With a 16-day revisit cycle, the Landsat images provide long, continuous and global records of the Earth's surface at a mesoscale (30 m × 30 m). In this study, we used atmospherically corrected Landsat 5/8 surface reflectance images covering our study area from 1984 to 2020, which are freely available on the Google Earth Engine (GEE) cloud platform [31]. For Landsat 5, atmospheric distortion is corrected using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm [32], while for Landsat 8, atmospheric effects are eliminated using the Land Surface Reflectance Code (LaSRC, [33]). These images consist of visible, near-infrared and shortwave-infrared bands processed to orthorectified surface reflectance and thermal infrared bands processed to orthorectified brightness temperature. For each image, cloud and cloud shadows were masked based on the CFMASK algorithm [34]. We excluded images in which none of the UGSs could be seen. Landsat 7 images were not used because of the data gaps caused by the failure of Figure 2. Overall workflow of the study. The main steps consist of data preprocessing, parameter retrieval and dynamic cooling assessment of UGSs. Note that SR refers to surface reflectance.

Landsat Surface Reflectance Collection
With a 16-day revisit cycle, the Landsat images provide long, continuous and global records of the Earth's surface at a mesoscale (30 m × 30 m). In this study, we used atmospherically corrected Landsat 5/8 surface reflectance images covering our study area from 1984 to 2020, which are freely available on the Google Earth Engine (GEE) cloud platform [31]. For Landsat 5, atmospheric distortion is corrected using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm [32], while for Landsat 8, atmospheric effects are eliminated using the Land Surface Reflectance Code (LaSRC, [33]). These images consist of visible, near-infrared and shortwave-infrared bands processed to orthorectified surface reflectance and thermal infrared bands processed to orthorectified brightness temperature. For each image, cloud and cloud shadows were masked based on the CFMASK algorithm [34]. We excluded images in which none of the UGSs could be seen. Landsat 7 images were not used because of the data gaps caused by the failure of the scan-line corrector. These processes finally gave us 287 Landsat 5  Directly combining Landsat 8 data with those from previous Landsat sensors is not suitable for long-term dynamic monitoring and trend analysis because Landsat 8 OLI/TIRS has spectrally narrower bands than those of previous Landsat sensors [35]. To improve the temporal continuity of the Landsat 8 data, we adjusted the surface reflectance of visible (blue, green, red), near-infrared and two shortwave bands of each Landsat 8 image to be compatible with Landsat 5 surface reflectance by using Roy's between-sensor transformation functions [36] where is the surface reflectance of Landsat 8 bands and is adjusted surface reflectance. The functions are derived using the ordinary least squares (OLS) regression from almost 60 million 30 m common pixel values extracted from more than 6000 Landsat 7 ETM+ and Landsat 8 OLI overlapped images acquired over three winter and three summer months for all the conterminous United States [36]. Other applications of the functions to long time series analyses of vegetation in the Tibetan Plateau [37] and Hindu Kush Himalaya [38] further indicate their universality.

Normalized Difference Vegetation Index
We employed the most widely used vegetation index, i.e., NDVI, to quantify vegetation greenness, density and conditions. NDVI is computed as the difference between the reflectance of the near-infrared and red bands, divided by their sum as follows: Directly combining Landsat 8 data with those from previous Landsat sensors is not suitable for long-term dynamic monitoring and trend analysis because Landsat 8 OLI/TIRS has spectrally narrower bands than those of previous Landsat sensors [35]. To improve the temporal continuity of the Landsat 8 data, we adjusted the surface reflectance of visible (blue, green, red), near-infrared and two shortwave bands of each Landsat 8 image to be compatible with Landsat 5 surface reflectance by using Roy's between-sensor transformation functions [36] as follows: Near in f rared 0.0306 + 0.8639ρ, Shortwave in f rared1 0.0116 + 0.9165ρ, Shortwave in f rared2 (1) where ρ is the surface reflectance of Landsat 8 bands and ρ adj is adjusted surface reflectance. The functions are derived using the ordinary least squares (OLS) regression from almost 60 million 30 m common pixel values extracted from more than 6000 Landsat 7 ETM+ and Landsat 8 OLI overlapped images acquired over three winter and three summer months for all the conterminous United States [36]. Other applications of the functions to long time series analyses of vegetation in the Tibetan Plateau [37] and Hindu Kush Himalaya [38] further indicate their universality.

Normalized Difference Vegetation Index
We employed the most widely used vegetation index, i.e., NDVI, to quantify vegetation greenness, density and conditions. NDVI is computed as the difference between the reflectance of the near-infrared and red bands, divided by their sum as follows:

Land Surface Temperature
We employed a generalized single-channel algorithm [39] to retrieve LST from atsensor brightness temperature using the following equations: where ε is the surface emissivity that is first spectrally adjusted from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global emissivity dataset and then modified for vegetation phenology and snow-cover changes following the procedures given by Malakar et al. [40]; L sen is the at-sensor radiance converted from raw Landsat digital numbers (DNs); T sen is the brightness temperature; λ is the effective wavelength; c 1 and c 2 are Planck's radiation constants; and ψ 1 , ψ 2 and ψ 3 are the atmospheric functions of total atmospheric water vapor content (w), given by the following: where c ij values are fitted using atmospheric sounding databases for different Landsat sensors [41]. The total atmospheric water vapor content was extracted from the Moderate Resolution Imaging Spectroradiometer (MODIS) combined multiangle implementation of atmospheric correction (MAIAC) land aerosol optical depth (AOD) gridded Level 2 product (MCD19A2) with a spatial resolution of 1 km [42]. Because the dataset is available after 2000, we used the water vapor from the National Centers for Environmental Prediction (NCEP) and National Center for Atmospheric Research (NCAR) reanalysis data as an alternative before 2000. When the MODIS water vapor was masked due to cloud contamination, the NCEP/NCAR data were used. The NCEP/NCAR water vapor has a 6-h temporal resolution (00:00, 06:00, 12:00 and 18:00 UTC) and 2.5-degree spatial resolution [43]. The water vapor at each Landsat acquisition time was approximated by linearly interpolating the two closest water vapor values [44].
We validated the retrieved LST with cloud cover less than 30% using in situ measurements of six surface radiation budget (SURFRAD) stations in the United States ( Figure 4). The stations primarily record the downwelling and upwelling components of broadband solar radiation and thermal infrared irradiance in 3 min intervals before January 2009 and every minute thereafter. The data are quality controlled and cover multiple land cover types (e.g., vegetation, grassland and cropland) in diverse climates [40]. The measurements used to validate LST were from 1995-2020. To minimize the uncertainty caused by spatial heterogeneity, we removed the Landsat LSTs that had a standard deviation of 3 × 3 pixels surrounding the station larger than 1 • C [40]. We also removed potential outliers based on the robust "3σ-Hampel identifier" [44][45][46]. These processes finally resulted in a total of

Evapotranspiration
We transferred the MODIS terrestrial ET algorithm [47] to retrieve fine-resolution ET by using the Landsat-based fractional vegetation cover. The algorithm is based on the Penman-Monteith equation [48]: where is the slope of the saturation water vapor pressure curve; is the available energy partitioned between sensible heat, latent heat and soil heat fluxes on the land surface; is the mean air density at constant pressure; is the specific heat of the air; is the saturation vapor pressure; is the actual water vapor pressure; − represents the vapor pressure deficit of the air; is the psychometric constant and and are the aerodynamic and surface resistance, respectively. The derivation of these inputs requires climate variables such as daily surface net radiation, air temperature and relative humidity, which were extracted from the NCEP climate forecast system version 2 (CFSv2) 6-h products [49].
Then, the total daily ET is the sum of the evaporation from the wet canopy surface, vegetation transpiration and evaporation from the soil surface calculated from Equation (7). When the relative humidity is less than 70%, the surface is assumed to be not covered by water, and the wet canopy surface evaporation and soil surface evaporation are then zero. A detailed calculation of these components can be found in [47].

Evapotranspiration
We transferred the MODIS terrestrial ET algorithm [47] to retrieve fine-resolution ET by using the Landsat-based fractional vegetation cover. The algorithm is based on the Penman-Monteith equation [48]: where s is the slope of the saturation water vapor pressure curve; A is the available energy partitioned between sensible heat, latent heat and soil heat fluxes on the land surface; ρ is the mean air density at constant pressure; c p is the specific heat of the air; e sat is the saturation vapor pressure; e is the actual water vapor pressure; e sat − e represents the vapor pressure deficit of the air; γ is the psychometric constant and r a and r s are the aerodynamic and surface resistance, respectively. The derivation of these inputs requires climate variables such as daily surface net radiation, air temperature and relative humidity, which were extracted from the NCEP climate forecast system version 2 (CFSv2) 6-h products [49]. Then, the total daily ET is the sum of the evaporation from the wet canopy surface, vegetation transpiration and evaporation from the soil surface calculated from Equation (7). When the relative humidity is less than 70%, the surface is assumed to be not covered by water, and the wet canopy surface evaporation and soil surface evaporation are then zero. A detailed calculation of these components can be found in [47].
In the MODIS algorithm, parameters that affect stomatal conductance, such as the mean potential stomatal conductance per unit leaf area, were parameterized globally using 46 AmeriFlux tower-based ET for 11 different biomes, i.e., evergreen needleleaf forest; evergreen broadleaf forest; deciduous needleleaf forest; deciduous broadleaf forest; mixed forest; woody savannas; savannas; closed shrubland; open shrubland; grassland, urban and built-up land; barren or sparsely vegetated land; and cropland. To determine which set of these parameters should be used for the UGSs in the study, we merged a high-resolution (100 m) land-cover scheme from the Copernicus global land-cover layers [50] in 2015 to correspond to the 11 biomes above. The land cover types of the permanent UGSs in this study are assumed to be unchanged during the study period.

Albedo
We converted narrowband surface reflectance to broadband albedo in total shortwave using the formula given by Liang [51] as follows: where α 1 , α 3 , α 4 , α 5 and α 7 are the narrowband surface reflectance of the Landsat blue band, red band, near-infrared band, shortwave-infrared band 1 and shortwave-infrared band 2, respectively. The MODIS daily albedo product (MCD43A3) was successfully used for change detection [52] and thus we employed it to help detect the albedo dynamics of UGSs.

Surface-Cooling Dynamics and Time Series Analysis
The surface-cooling effects of UGSs at time t are defined as the LST difference between UGSs and Tian'anmen Square (Figure 1), namely, Tian'anmen Square was chosen as a typical urban impervious surface for comparison because of its homogeneity and invariance since 1954. Given this, a negative ∆LST indicates that the surface temperature of UGSs is lower than that of urban impervious surfaces; therefore, they represent cooling effects and smaller values suggest stronger cooling.
To detect and explain the dynamic cooling effects of UGSs, we first generated four time series for the NDVI, ∆LST, ET and albedo of the selected permanent UGSs (Figure 1) from 1984-2020. Then, we imputed the missing data (due to the removal of cloudy images) of each time series by using the seasonal Kalman filter and deconstructed the four completed time series into an additive seasonal component (S t ), trend component (T t ) and residual component (E t ) using a locally weighted regression smoother (LOESS) [53], as follows: where the seasonal component represents the vegetation phenology cycle, and the trend component detects the temporal vegetation dynamics. Hereafter, we implemented an iterative procedure of linear fitting [54] on the four trend series to examine potential breakpoints and whether segmented relationships existed in the time series.

Greenness Dynamics
Based on the reconstructed Landsat NDVI time series from 1984 to 2020 (Figure 5b), we extracted the NDVI trend of UGSs and observed overall greening (Figure 5d). In 1986, even in the growing seasons, UGSs were in poor condition with a very low NDVI of approximately 0.2 (Figure 5b), and they were hardly distinguished from non-vegetated areas (Figure 5f1). By 2020, however, the annual UGS NDVI had risen to 0.4 (Figure 5d), with the growing season NDVI reaching 0.5 (Figure 5b). Such changes were achieved through three segmented relationships before 1996, during 1996-2000 and after 2000. Before 1996, the NDVI slowly increased at a rate of +0.004/year After 1996, however, the UGSs experienced short-term deforestation (with an NDVI decrease rate of −0.012/year) until 2000, reducing their greenness almost back to the level of 1984. Hereafter, UGSs rapidly grew with a sharply increased NDVI rate of +0.013/year, which was particularly true for large UGSs (circled in Figure 5f).
Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 24 with a sharply increased NDVI rate of +0.013/year, which was particularly true for large UGSs (circled in Figure 5f).

Responses of Urban Green Space Surface Cooling to Greening
By comparing the temporal LST difference (∆ ) between UGSs and Tian'anmen Square (a typical permanent impervious surface in Beijing, Figure 1), we observed a positive response of most UGS daytime surface-cooling effects (∆ ) to their greenness (Figures 5 and 6). UGS greening (before 1996 and after 2000) led to cooling enhancements, while their browning (during 1996-2000, Figure 5d) resulted in cooling attenuation (Figure 6c). The overall cooling enhancement made UGSs recently experience an approximately 2 °C lower annual surface temperature than impervious surfaces compared to the value of no more than 1.5 °C in 1984 (Figure 6c). A spatial comparison of this cooling enhancement was illustrated by using the ∆ image composites from 2000 and 2020. This is because 2000 was a trough representing poor vegetation performance, which could be sharply compared with 2020, when the vegetation condition had grown to its peak. In contrast with 1984, 2000 also had sufficient images to make an annual composite ( Figure  6a). Spatially, almost all UGSs experienced a distinctly lower LST in 2020 (Figure 6e), which was true only for approximately half of the UGSs in 2000 (Figure 6d).

Responses of Urban Green Space Surface Cooling to Greening
By comparing the temporal LST difference (∆LST) between UGSs and Tian'anmen Square (a typical permanent impervious surface in Beijing, Figure 1), we observed a positive response of most UGS daytime surface-cooling effects (∆LST) to their greenness ( Figures 5 and 6). UGS greening (before 1996 and after 2000) led to cooling enhancements, while their browning (during 1996-2000, Figure 5d) resulted in cooling attenuation (Figure 6c). The overall cooling enhancement made UGSs recently experience an approximately 2 • C lower annual surface temperature than impervious surfaces compared to the value of no more than 1.5 • C in 1984 (Figure 6c). A spatial comparison of this cooling enhancement was illustrated by using the ∆LST image composites from 2000 and 2020. This is because 2000 was a trough representing poor vegetation performance, which could be sharply compared with 2020, when the vegetation condition had grown to its peak. In contrast with 1984, 2000 also had sufficient images to make an annual composite (Figure 6a). Spatially, almost all UGSs experienced a distinctly lower LST in 2020 (Figure 6e), which was true only for approximately half of the UGSs in 2000 (Figure 6d).
x FOR PEER REVIEW 10 of 24 We also found that UGS daytime surface cooling significantly increased (p < 0.05, ttest) to varying degrees in different seasons except for winter (Figure 7). The most noticeable cooling enhancement occurred on summer days (Figure 7a2-c2). Although vegetation provides strong cooling effects in summer, UGSs in the summer of 2000 did not show much lower ∆ values than those in spring and autumn, with a mean ∆ of approximately −2 °C. However, this was not the case in the summer of 2020, when the distribution of ∆ integrally shifted to the left (Figure 7c2), and the mean cooling effects increased to more than 4 °C. This change was also illustrated by the fact that 88% of USG pixels represented a negative ∆ in 2000, while the proportion had increased to 97% in 2020 (Figure 7c2). For spring and autumn, although the boundaries of ∆ remained almost unchanged during the two periods (Figure 7c1,c3), we observed evident left movements of the wave crests that resulted in mean cooling effect increases of 0.7 °C for spring and 0.5 °C for autumn. In winter, we failed to observe a significant cooling difference for UGSs. UGSs in both 2000 and 2020, on average, had a near-zero positive ∆ ( Figure  7c4), implying their dim warming effects. We also found that UGS daytime surface cooling significantly increased (p < 0.05, t-test) to varying degrees in different seasons except for winter (Figure 7). The most noticeable cooling enhancement occurred on summer days (Figure 7a2-c2). Although vegetation provides strong cooling effects in summer, UGSs in the summer of 2000 did not show much lower ∆LST values than those in spring and autumn, with a mean ∆LST of approximately −2 • C. However, this was not the case in the summer of 2020, when the distribution of ∆LST integrally shifted to the left (Figure 7c2), and the mean cooling effects increased to more than 4 • C. This change was also illustrated by the fact that 88% of USG pixels represented a negative ∆LST in 2000, while the proportion had increased to 97% in 2020 ( Figure 7c2). For spring and autumn, although the boundaries of ∆LST remained almost unchanged during the two periods (Figure 7c1,c3), we observed evident left movements of the wave crests that resulted in mean cooling effect increases of 0.7 • C for spring and 0.5 • C for autumn. In winter, we failed to observe a significant cooling difference for UGSs. UGSs in both 2000 and 2020, on average, had a near-zero positive ∆LST (Figure 7c4), implying their dim warming effects.

Biophysical Mechanisms
To explain the change in UGS surface-cooling effects, we established the dynamics of two main biophysical parameters (i.e., ET and albedo) that largely regulated vegetation heat exchange. ET characterizes the capability of vegetation to transfer sensible heat into latent heat, and vegetation with higher ET provides more cooling effects (transpiration cooling). In contrast, albedo represents how much incidental radiation can be absorbed by vegetation, and low albedo values imply that more energy will be absorbed to produce warming effects. Thus, a combined effect of transpiration cooling and albedo warming determines the final thermal responses to changes in vegetation conditions. As UGSs continued to turn greener, their ET and albedo values were positively and negatively responsive, respectively (

Biophysical Mechanisms
To explain the change in UGS surface-cooling effects, we established the dynamics of two main biophysical parameters (i.e., ET and albedo) that largely regulated vegetation heat exchange. ET characterizes the capability of vegetation to transfer sensible heat into latent heat, and vegetation with higher ET provides more cooling effects (transpiration cooling). In contrast, albedo represents how much incidental radiation can be absorbed by vegetation, and low albedo values imply that more energy will be absorbed to produce warming effects. Thus, a combined effect of transpiration cooling and albedo warming determines the final thermal responses to changes in vegetation conditions. As UGSs continued to turn greener, their ET and albedo values were positively and negatively responsive, respectively ( Such competing biophysical changes resulted in segmented changes in UGS surface cooling (Figure 8b2). Due to the opposite thermal effects of ET and albedo, transpiration cooling dominated albedo warming in both the two greening periods (1984-1996 and 2000-2020) to provide discernable cooling enhancements. However, we found that a higher greening rate resulted in a lower cooling increase rate since the increase rate of UGS surface cooling after 2000 (−0.018 °C/year, Figure 8b2) was even less than half of that before 1996 (−0.050 °C/year, Figure 8b2). This could be related to the fact that albedo is more sensitive to UGS greening than ET with a comparison of the relative decrease rate of albedo (( − )/ = 85%, Figure 8d2) and the relative increase rate of ET (( − )/ = 8%, Figure 8c2) in the two greening stages. Therefore, the cooling increase rate of UGS in the latter stage was largely offset by albedo-induced warming.

Greening of Urban Green Spaces
The highlighted dynamics of UGS daytime cooling for the long run closely correlate with UGS vegetation conditions. Our results showed an overall increase in permanent UGS greenness during the past 40 years (1984-2020, Figure 5), implying that vegetation, at least in large UGSs, is healthier and more ecologically functional than it was before. This result may not be surprising because mainly due to CO2 fertilization, nitrogen deposition and climate changes, most of the planet has observed a phenomenon of 'greening' in recent decades [55,56]. High-intensity urbanization leads to higher radiation forces and temperatures, further enhancing vegetation growth and prolonging vegetation phenology to advance/delay the start/end of the growing season [57,58]. Furthermore, most Chinese UGSs are managed by municipal governments, and thus, they are more frequently influenced by human activities than UGSs in western cities [30]. After 2000, the Beijing Municipal Forestry and Parks Bureau implemented extensive park improvements, such as planting/replanting, transplanting, maintenance, protection from fire and diseases and the introduction of nonnative species. The activities play a critical role in vegetation (e.g., trees, shrubs, herbaceous plants) conditions and overall greenery coverage. For instance, more sophisticated management practices, e.g., better irrigation increases soil moisture to trigger vegetation greening [24], which is especially true for water-limited drylands. Other adverse effects, such as particulate matter accumulation, may decrease the photosynthetic Such competing biophysical changes resulted in segmented changes in UGS surface cooling (Figure 8b2). Due to the opposite thermal effects of ET and albedo, transpiration cooling dominated albedo warming in both the two greening periods (1984-1996 and 2000-2020) to provide discernable cooling enhancements. However, we found that a higher greening rate resulted in a lower cooling increase rate since the increase rate of UGS surface cooling after 2000 (−0.018 • C/year, Figure 8b2) was even less than half of that before 1996 (−0.050 • C/year, Figure 8b2). This could be related to the fact that albedo is more sensitive to UGS greening than ET with a comparison of the relative decrease rate of albedo ((k 2 − k 1 )/k 1 = 85%, Figure 8d2) and the relative increase rate of ET ((k 2 − k 1 )/k 1 = 8%, Figure 8c2) in the two greening stages. Therefore, the cooling increase rate of UGS in the latter stage was largely offset by albedo-induced warming.

Greening of Urban Green Spaces
The highlighted dynamics of UGS daytime cooling for the long run closely correlate with UGS vegetation conditions. Our results showed an overall increase in permanent UGS greenness during the past 40 years (1984-2020, Figure 5), implying that vegetation, at least in large UGSs, is healthier and more ecologically functional than it was before. This result may not be surprising because mainly due to CO 2 fertilization, nitrogen deposition and climate changes, most of the planet has observed a phenomenon of 'greening' in recent decades [55,56]. High-intensity urbanization leads to higher radiation forces and temperatures, further enhancing vegetation growth and prolonging vegetation phenology to advance/delay the start/end of the growing season [57,58]. Furthermore, most Chinese UGSs are managed by municipal governments, and thus, they are more frequently influenced by human activities than UGSs in western cities [30]. After 2000, the Beijing Municipal Forestry and Parks Bureau implemented extensive park improvements, such as planting/replanting, transplanting, maintenance, protection from fire and diseases and the introduction of nonnative species. The activities play a critical role in vegetation (e.g., trees, shrubs, herbaceous plants) conditions and overall greenery coverage. For instance, more sophisticated management practices, e.g., better irrigation increases soil moisture to trigger vegetation greening [24], which is especially true for water-limited drylands. Other adverse effects, such as particulate matter accumulation, may decrease the photosynthetic rate to cause the poor performance of UGSs [59]. Indeed, a short-term greenness decline in 1996-2000 ( Figure 5d) revealed a trough in UGSs, which was consistent with the contemporaneous field surveys conducted in the central part of Beijing that indicate trees in 2002 were mainly dominated by younger age classes and were in poor condition [60].

The Necessity of Dynamic Monitoring of Urban Green Space Cooling Effects
Because many of these possible driving forces vary with urban development and climate change, long-term monitoring becomes particularly important for urban vegetation eco-functional assessments. However, most existing conclusions related to UGS cooling effects have been drawn from contrast experiments conducted over a few years [11,[61][62][63] or even on a single day [64,65]. Short-term observations might cause an over-or underestimation of long-term cooling potential to varying degrees, depending on how fast the urban area grows and how well UGSs are managed, as well as their interactions with changing environments. In this study, UGSs in the past 40 years suggested surface-cooling enhancements in all seasons ranging from less than 1 • C to more than 2 • C. On warm summer days when UGSs are most needed to provide cooling, the large enhancements not only further reduce the surface temperature, but contribute to the formation of a convection cycle to provide canopy-layer cooling by subsiding the warmer urban air from above into UGSs and blowing their cooler air [66].
The cooling increase rate of UGS is not as fast as its greening rate because increasing transpiration cooling was largely offset by albedo-induced energy absorption, given the sharp decline of albedo since 2000 (Figure 8d2). This implies that more irrigation water might have been used to fulfill the increasing canopy transpiration demands, while the annually averaged surface-cooling benefits compared to 1984 were observed to be less than 1 • C (Figure 6c). Although this study was conducted on permanent UGSs, the results have important implications for urban afforestation, since these changes in urban vegetation physiology could be harbingers of newly built green belts and urban parks and those areas that will be a large area of green spaces in the next decade. According to the Beijing urban master plan (2016-2035), multiple urban afforestation and greening programs will be conducted to increase per capita green space [12], and the city's overall forest cover rate should increase to no less than 45% in 2035. Therefore, urban planners and policymakers should consider the dynamics of the biophysical features of vegetation found in this study to balance the maintenance costs and the environmental gains in future UGS management.

Limitations and Future Work
The use of a series of Landsat-based datasets in this study is mainly due to its long temporal coverage (since 1984) and capacity to discern human-scale processes, such as the managed UGSs in Beijing. Because of its instantaneous imaging, which occurs at the satellite overpass time (~10:53 a.m. local time in Beijing), our results represent the long-term dynamics of the biophysical characteristics of vegetation close to midday when illumination and water vapor reach their maximum and minimum, respectively. Thus, they inevitably sacrifice diurnal changes in ∆LST and thus, UGS cooling effects. Additionally, due to the lack of nighttime high-resolution thermal infrared data, the temporal changes in ∆LST are not evaluated at night when the urban heat island (UHI) intensity is generally strong. In addition, both LST and air temperature have been commonly used as key indicators measuring the UHI intensity or UGS cooling effects. However, despite high correlation, they are hardly comparable at the city scale because LST is the radiative skin temperature of the land surface that mainly depends on land use and land cover. In contrast, air temperature represents the thermal conditions of urban canopy layers, which are largely affected by geometry-induced advection [67]. Accordingly, future studies focusing on the dynamics of air temperature differences inside and outside growing UGSs are highly encouraged.

Conclusions
Long-term dynamic monitoring plays a key role in the development and management of UGSs. Revealed by Landsat time series data on the past 40 years, we found segmented changes in UGS cooling effects that were mainly linked to the responses of canopy transpiration and albedo to vegetation conditions. During a rapid greening of UGSs in the recent two decades, transpiration cooling dominated albedo-induced warming to provide a discernable cooling enhancement, despite relatively more sensible changes in albedo. Such enhancement showed seasonal differences ranging from less than 1 • C to more than 2 • C, and the most evident enhancements occurred on summer days (~2.4 • C). These findings extend the current understanding of UGS cooling effects and provide insights into how changes in vegetation physiology affect urban climate by altering the surface energy balance and the hydrologic cycle.

Acknowledgments:
We thank the three anonymous reviewers, academic editor and Qiang Liu at Beijing Normal University for their valuable comments on the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.