**1. Introduction**

Owing to the high heterogeneity and complexity in urban ecosystems, it is rather difficult to monitor or predict the hydrological dynamics of urban surfaces [1]. Some megacities, e.g., Beijing—the capital city of China—have experienced strong urbanization, large population inflow, island effect and climate change during the past few decades [2]. These changes induce urban hydrological processes to be highly uncertain and make policymakers face tough challenges in water use planning and management. Therefore, there is an urgen<sup>t</sup> need to accurately estimate urban hydrological processes.

Evapotranspiration (ET), as a key component of the urban hydrological processes and surface energy balance, plays an important role in regulating water resource supply and relieving the urban island effect (e.g., surface cooling) [3]. Different from natural ecosystems, the urban ecosystems include large proportions of artificial modifications in land cover, such as impervious surfaces including roofs, squares and cement or asphalt roads. These man-made reconstructions could contribute a large fraction of evaporation [4–6], but the quantification at city levels remains highly uncertain due to a lack of clearly distinguishing estimations of ET between impervious surfaces and vegetated or bare-soil lands. The

*Article*

**Citation:** Zhang, X.; Song, P. Estimating Urban Evapotranspiration at 10m Resolution Using Vegetation Information from Sentinel-2: A Case Study for the Beijing Sponge City. *Remote Sens.* **2021**, *13*, 2048. https:// doi.org/10.3390/rs13112048

Academic Editor: José A. Sobrino

Received: 14 May 2021 Accepted: 21 May 2021 Published: 22 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

good news is that recent developments of fine resolution remote sensing for land use and land cover classification, vegetation dynamics and environmental monitoring provide new opportunities to estimate urban ET more accurately [7]. For example, the Sentinel-2, as part of the European Commission's Copernicus program with the launch of satellite Sentinel-2A on 23 June 2015, are monitoring variability in global land surface conditions at a 10–60m resolution and a 5–10-day revisit [8].

In this study, we take the sponge city project in Beijing city as a study case to estimate ET of the urban ecosystem at 10m resolution using the satellite-based land cover map and vegetation information derived from Sentinel-2 data. Beijing city has 5-year mean annual precipitation of 560 mm and mean annual temperature of 12 ◦C, with the potential evaporation of 550~600 mm year<sup>−</sup>1. To reduce stress on water supply (e.g., ~30 m<sup>3</sup> per capita water use per year in Beijing) and urban environment, the Beijing Sponge City project was started on 4 December 2017, aiming at turning 20% of Beijing city into a sponge city covered area by 2020 [9]. Therefore, to evaluate the benefit of this project, it is essential to implement a city-level assessment of the project-induced ecohydrological changes at fine resolution.

#### **2. Materials and Methods**

#### *2.1. Observational Forcing Datasets*

2.1.1. Land Cover Map at 10m Resolution Derived from Sentinel-2

The land cover classification (LCC) global map at 10m resolution was obtained from FROM-GLC10 [7]. The FROM-GLC10 LCC data is developed based on Sentinel-2 data in 2017 with Google Earth Engine, and the overall accuracy of this LCC validated against the circa 2015 validation sample is 73% [7]. The LCC data includes 10 classes (i.e., cropland, forest, grassland, shrubland, wetland, tundra, impervious surface, bare land, and snow/ice). The most advances of the FROM-GLC10 LCC map compared to previous Landsat series-based LCC products are that it provides more spatial detail, better distinguish the forest from shrub or grassland classes, and better performance in coastal areas [7].

#### 2.1.2. NDVI and LAI at 10m Resolution Derived from Sentinel-2

We calculate the Normalized Difference Vegetation Index (NDVI) from the Sentinel-2 reflected radiance by

$$NDVI = \frac{R\_{nir} - R\_{rad}}{R\_{nir} + R\_{red}} \tag{1}$$

where *Rnir* and *Rred* are the spectral bands at near infrared (842 nm) and red (665 nm), respectively. The Sentinel-2 reflectance data are available at the USGS EROS Center (https://www.usgs.gov/centers/eros, accessed on 8 April 2021). The leaf area index (LAI) at a 10m resolution was derived from the retrieved Sentinel-2 NDVI using a nonlinear regression model between LAI and NDVI,

$$LAI = a \ast \exp(b \ast NDVI) + c \tag{2}$$

where the parameters *a*, *b* and *c* are determined as 12.4, 6.4 and 0.6, respectively.

The determination process was based on MODIS-based NDVI and LAI products, which was described as: (i) The MODIS LAI (MOD15A2H) and surface reflectance (SR) products (MOD09A1) at 500m were collected over the study area (Beijing Sponge City) for the summer months (June, July, and August) from 2013 to 2019. The NDVI was then estimated using the 500m SR product (Equation (1)). (ii) MODIS LAI and NDVI values were collocated on the pixel basis. As the MODIS LAI product has a valid range between 0 and 6.9, with a precision of 0.1, we classified all NDVI values into 69 groups based on unique LAI values (eliminating the zero-LAI group). The probability density plots for each group are shown in Figure 1a. (iii) For each LAI-based value group, the probability density was fitted using the Gaussian distribution function. Then, the NDVI value corresponded by the maximum probability density was extracted and collocated with the specific LAI value.

The result shows a strong exponential relationship between NDVI and LAI, especially when the LAI value increased beyond 0.5–0.6 (Figure 1b). Therefore, the scatter values were fitted using the exponential model (Equation (2)), which resulted in an R<sup>2</sup> of 0.82, implying that such exponential model in Equation (2) is robust for the Beijing Sponge City.

**Figure 1.** The nonlinear relationship between LAI and NDVI. (**a**) Frequency distribution of NDVI–LAI. Higher values of the scatter number in (**a**) indicate stronger relation between LAI and NDVI. (**b**) Regression between LAI and NDVI. The brown curve was fitted using Equation (2) from the NDVI–LAI values (blue stars) under the maximum probability density.

#### 2.1.3. Surface Climate Driving Dataset

To estimate evapotranspiration of urban ecosystem at 10m resolution, a high-resolution surface climate forcing data including precipitation, surface air temperature, wind speed, surface pressure, specific humidity, downward shortwave and longwave radiations, etc., is needed to drive the terrestrial evapotranspiration model (PML-V2 model, see Section 2.2). In this study, we used the China Meteorological Forcing Dataset (CMFD) version 1 at 0.1◦ × 0.1◦ and daily resolution for June 2018 as input for the PML-V2 model. The CMFD V1.0 dataset covered the period of 1979–2018 and was downscaled from station-based data, TRMM satellite-based precipitation, GEWEX-SRB shortwave radiation and the GLDAS forcing dataset [10]. The surface climate driving variables used for the Beijing Sponge City area were spatially bilinearly interpolated onto a 10m × 10m resolution. The monthly CO2 concentration observed in June 2018 is set as 407 ppm.

## *2.2. PML-V2.1 Model*

Version 2 of the Penman–Monteith–Leuning model (PML-V2) was developed by coupling the widely-used photosynthesis model [11] and a canopy stomatal conductance model [12] with the Penman–Monteith energy balance equation [13] to jointly estimate gross primary productivity (GPP), *Ec* and *Es* [14–18]. The PML-V2 model also simulates the *Ei* based on a revised Gash-model scheme [19]. The PML-V2 model has been applied to successfully produce the MODIS LAI-based global GPP and ET products at a 500m and 8-day resolution from 2002 to present, which were noticeably better than most widely used GPP and ET products [16]. In this study, we incorporated modules of impervious surface evaporation (*Eu*) and open-water evaporation (*Ew*) into the PML model (PML-V2.1) to make it suitable for urban ecosystems. Key parameters used in the PML-V2.1 model are provided in Table 1. The following shows the detailed description for PML-V2.1.


**Table 1.** Key parameters used in the PML-V2.1 model.

(a) CRO: cropland, MIF: mixed forest, GRA: grassland, SHR: shrubland, WET: wetland, WAT: water body, IMP: impervious surface and BAR: bare land.

#### 2.2.1. Energy Balance at Urban Land Surface

Based on the surface energy balance, the net radiation (*Rn*) can be balanced by the latent heat flux (*LE*), sensible heat flux (*H*) and ground heat flux (*G*). As for a biweekly or longer estimation, the *G* is often negligible (*G H* + *LE*), then the *H* is given by

$$H = R\_{\rm II} - LE - G \approx R\_{\rm II} - LE \tag{3}$$

The net radiation at the surface is the sum of the net shortwave downward radiation and the net longwave downward radiation,

$$R\_n = (1 - \alpha)SW + \left(LW - \epsilon \sigma T\_a^4 \right) \tag{4}$$

where the shortwave downward radiation *SW* (W m<sup>−</sup>2), the longwave downward radiation *LW* (W m<sup>−</sup>2) and the surface air temperature *Ta* (K) are from atmospheric forcing input data [10]. The shortwave albedo *α* (-) and the longwave emissivity (-) are from satellitebased estimations. *σ* is the Stefan–Boltzmann constant (5.67 × 10−<sup>8</sup> W m<sup>−</sup><sup>2</sup> <sup>K</sup>−4). The latent heat flux (*LE*,Wm−2) is calculated by *LE* = 1*c λET* and

$$ET = E\_c + E\_s + E\_i + E\_u + E\_w \tag{5}$$

where *λ* = 2500 − 2.2(*Ta* − 273.15) is the latent heat of vaporization (kJ kg−1) at *Ta*, and *c* (=86.4) is a conversion factor for units from (MJ m<sup>−</sup><sup>2</sup> <sup>d</sup>−1) to (W m<sup>−</sup>2). *ET* is the evapotranspiration (mm <sup>d</sup>−1) summed from the canopy transpiration (*Ec*) and soil evaporation (*Es*), interception evaporation (*Ei*), impervious surface evaporation (*Eu*) and open-water evaporation (*Ew*).

#### 2.2.2. Canopy Transpiration (*Ec* ) and Soil Evaporation (*E*s)

The transpiration at canopy scales (*Ec*) is coupled with the photosynthesis process (*Ags*) via the dynamical modulation of the canopy stomatal conductance (*Gc*), and the soil evaporation (*Es*) depends on absorbed energy flux and soil water deficit,

$$E\_{\mathfrak{c}} = \frac{\varepsilon A\_{\mathfrak{c}} + \left(\rho c\_p / \gamma\right) D\mathcal{G}\_{\mathfrak{u}}}{\varepsilon + 1 + \mathcal{G}\_{\mathfrak{u}} / \mathcal{G}\_{\mathfrak{c}}} \tag{6}$$

$$E\_s = \frac{f\varepsilon A\_s}{\varepsilon + 1} \tag{7}$$

where the surface available energy (*A* = *Rn* − *G*) is divided into canopy absorbed energy (*Ac*) and soil absorbed energy (*As*), *Ac* = (1 − *τ*)*<sup>A</sup>* and *As* = *<sup>τ</sup>A*, *τ* = *exp*(−*kALAI*), *kA* = 0.6. *ε* = Δ*γ* , and Δ is the slope of the curve relating saturation water vapor pressure to temperature (kPa ◦C−1). *ρ* is the density of air (g m<sup>−</sup>3); *cp* is the specific heat of air at constant pressure (J g<sup>−</sup><sup>1</sup> ◦C−1); *D* is vapor pressure deficit (kPa); *Ga* is the aerodynamic conductance (m s<sup>−</sup>1); *Gc* (m s<sup>−</sup>1) is the canopy conductance; *f* is a dimensionless variable that determines the water availability for soil evaporation.

The canopy stomatal conductance (*Gc*) is calculated by the photosynthesis process (*Ags*) for each PFT with the constraint of *D* at surface.

$$G\_{\mathfrak{c}} = \frac{mA\_{\mathfrak{g}^{\mathfrak{s}}}}{\mathbb{C}\_{\mathfrak{a}}(1 + D/D\_0)}\tag{8}$$

$$A\_{\mathcal{S}^s} = \frac{P\_1 \mathbb{C}\_{\mathfrak{u}}}{k(P\_2 + P\_4)} \left\{ kLAI + ln \frac{P\_2 + P\_3 + P\_4}{P\_2 + P\_3 \exp(kLAI) + P\_4} \right\} \tag{9}$$

$$V\_{\rm m} = \frac{V\_{m,25} \exp[a(T - 25)]}{1 + \exp[b(T - 41)]} \tag{10}$$

where *P*1 = *AmβI*0*η*; *P*2 = *AmβI*0; *P*3 = *AmηCa*; *P*4 = *βI*0*ηCa*; *Am* = 0.5*Vm*. *I*0 is the photosynthetically active radiation (PAR, in mol) from shortwave downward radiation. *Ca* is the atmospheric CO2 concentration (in ppm or mol mol−1). *Vm*,25, *β*, *η*, *m*, *Dmin*, *Dmax* and *D*0 are the parameters (see Table 1) in the PML-V2.1 model.

Finally, the gross primary productivity (*GPP*) is calculated by

$$GPP = A\_{\mathbb{R}^s} f\_D \tag{11}$$

$$f\_D = \begin{cases} 1, & D < D\_{\min} \\ \frac{D\_{\max} - D}{D\_{\max} - D\_{\min}}, & D\_{\min} < D < D\_{\max} \\ 0, & D > D\_{\max} \end{cases} \tag{12}$$

where *fD* is the *D* constraint function; *Dmin* and *Dmax* are the parameters (see Table 1).

#### 2.2.3. Interception Evaporation (*Ei*) by Canopy Vegetation

The rainfall interception evaporation (*Ei*) is calculated by the van Dijk model, which was developed by van Dijk and Bruijnzeel [19,20], who modified it from the original Gash model [21,22]. The modified Van Dijk model used in this study is described by

$$E\_i = \begin{cases} \begin{array}{c} f\_V P \\ f\_V P\_{\text{wet}} + f\_{\text{ER}}(P - P\_{\text{wet}}) \end{array} & \text{for } P < P\_{\text{wet}} \\\end{array} \tag{13}$$

with

$$f\_V = 1 - \exp\left(-\frac{LAI}{LAI\_{ref}}\right) \tag{14}$$

$$P\_{\text{wet}} = -\ln\left(1 - \frac{f\_{\text{ER}}}{f\_V}\right)\frac{S\_V}{f\_{\text{ER}}}\tag{15}$$

$$f\_{ER} = f\_V F\_{ER0} \tag{16}$$

$$S\_V = S\_{leaf}LAI \tag{17}$$

where *P* is rainfall rate (mm <sup>d</sup>−1), and *Pwet* is reference threshold precipitation rate when the canopy is wet (mm <sup>d</sup>−1). *fV* describes the fraction area covered by intercepting leaves, which is determined by the leaf area index (*LAI*) and reference *LAI* (*LAIref* , see Table 1) for the specific vegetation types. *fER* is the ratio of the average evaporation rate over the average rainfall intensity during storms, and *SV* is the canopy rainfall storage capacity (mm), which currently is parameterized as water storage in the leaf at the canopy level. The *fER*0 and *Sleaf* are the parameters shown in Table 1.

#### 2.2.4. Impervious Surface Evaporation (*Eu*)

Impervious surface evaporation (*Eu*) is calculated by the revised van Dijk model in this study,

$$E\_{\rm II} = \begin{cases} \begin{array}{c} f\_{\rm II}P \\ f\_{\rm II}P\_{\rm wet} + f\_{\rm ER}(P - P\_{\rm wet}) \end{array} & \text{for } P < P\_{\rm wet} \\\end{cases} \tag{18}$$

with

$$f\_{l\bar{l}} = 1 - f\_V - f\_w \tag{19}$$

$$P\_{\text{wet}} = -\ln\left(1 - \frac{f\_{ER}}{f\_{U}}\right) \frac{S\_{U}}{f\_{ER}}\tag{20}$$

$$f\_{ER} = f\_{ll} F\_{ER0} \tag{21}$$

where *fU* describes the fractional area covered by impervious surface in urban ecosystems, which is the rest fraction of vegetation coverage (*fV*) and water body covered fraction (*fw*). *SU* is the impervious surface storage capacity (mm).

2.2.5. Open-Water Evaporation (*Ew*)

The open-water evaporation (*Ew*) in lakes, rivers and other water bodies for Northern China is parameterized on the basis of Dalton's Law [23],

$$E\_w = 0.144(1 + 0.75 \mathcal{U} I\_{1.5}) \left[ D + \Delta(T\_{1.5})(a - 1)T\_{1.5} \right] \tag{22}$$

where *U*1.5 and *T*1.5 is wind speed (m/s) and air temperature (◦C) at 1.5 m height, respectively. <sup>Δ</sup>(*<sup>T</sup>*1.5) is the slope of the curve relating saturation water vapor pressure to temperature *T*1.5. *α* − 1 is a regulator coefficient for atmospheric stability.
