Articles | Volume 16, issue 15
Model evaluation paper
10 Aug 2023
Model evaluation paper |  | 10 Aug 2023

Simulating heat and CO2 fluxes in Beijing using SUEWS V2020b: sensitivity to vegetation phenology and maximum conductance

Yingqi Zheng, Minttu Havu, Huizhi Liu, Xueling Cheng, Yifan Wen, Hei Shing Lee, Joyson Ahongshangbam, and Leena Järvi

The Surface Urban Energy and Water Balance Scheme (SUEWS) has recently been introduced to include a bottom-up approach to modeling carbon dioxide (CO2) emissions and uptake in urban areas. In this study, SUEWS is evaluated against the measured eddy covariance (EC) turbulent fluxes of sensible heat (QH), latent heat (QE), and CO2 (FC) in a densely built neighborhood in Beijing. The model sensitivity to maximum conductance (gmax) and leaf area index (LAI) is examined. Site-specific gmax is obtained from observations over local vegetation species, and LAI parameters are extracted by optimization with remotely sensed LAI obtained from a Landsat 7 data product. For the simulation of anthropogenic CO2 components, local traffic and population data are collected. In the model evaluation, the mismatch between the measurement source area and simulation domain is also considered.

Using the optimized gmax and LAI, the modeling of heat fluxes is noticeably improved, showing higher correlation with observations, lower bias, and more realistic seasonal dynamics of QE and QH. The effect of the gmax adjustment is more significant than the LAI adjustment. Compared to heat fluxes, the FC module shows lower sensitivity to the choices of gmax and LAI. This can be explained by the low relative contribution of vegetation to the net FC in the modeled area. SUEWS successfully reproduces the average diurnal cycle of FC and annual cumulative sums. Depending on the size of the simulation domain, the modeled annual accumulated FC ranges from 7.4 to 8.7 kgCm-2yr-1, compared to 7.5 kgCm-2yr-1 observed by EC. Traffic is the dominant CO2 source, contributing 59 %–70 % to the annual total CO2 emissions, followed by human metabolism (14 %–18 %), buildings (11 %–14 %), and CO2 release by vegetation and soil respiration (6 %–10 %). Vegetation photosynthesis offsets only 5 %–10 % of the total CO2 emissions. We highlight the importance of choosing the optimal LAI parameters and gmax when SUEWS is used to model surface fluxes. The FC module of SUEWS is a promising tool in quantifying urban CO2 emissions at the local scale and therefore assisting in mitigating urban CO2 emissions.

1 Introduction

Currently, half of the global population resides in urban areas, and this percentage is projected to grow to 68 % by the middle of the 21st century (United Nations Department of Economic and Social Affairs2019). Urban expansion has reshaped the morphological, thermal, and dynamical properties of the land surface (Grimmond and Oke2006; Oke1995; Zhu et al.2016). In addition, intensive human activities in urban areas have caused a large quantity of greenhouse gas emissions (Marcotullio et al.2013; Velasco and Roth2010). Both factors have influenced urban climate from micro-scales to regional scales (Johansson and Emmanuel2006; Sarangi et al.2018; Tan et al.2010). Climatic and environmental risks due to urbanization are frequently reported, such as heat waves, flooding, and air pollution (Qian et al.2022; Watts et al.2015). In this context, there is a pressing need to better understand the effects of urbanization on land–atmosphere interaction, preferably in a quantitative way.

Urban land surface models (ULSMs) are widely used to simulate urban land–atmosphere interactions, including the exchanges of energy, water, and CO2 and hydrological processes (Chen et al.2011; Masson et al.2013). The results from the First International Urban Land Surface Model Comparison Project suggested that the most important processes for urban surface energy balance were radiative and vegetation processes (e.g., vegetation fraction, seasonal cycle of vegetation phenology) (Best and Grimmond2015; Grimmond et al.2010; Nordbo et al.2015). Long-term observations with low vegetation cover (< 30 %) were especially needed to evaluate heat flux simulation, as energy distribution was found to be sensitive in such environments (Best and Grimmond2016).

The Surface Urban Energy and Water balance Scheme (SUEWS) is one of the widely tested ULSMs (Järvi et al.2011, 2014; Ward et al.2016). SUEWS is designed to run with surface information (e.g., surface cover fractions) and a minimum amount of model forcing data. In recent years, Supy (SUEWS in Python) was developed to allow Python front-end implementation for broader and easier applications (Sun and Grimmond2019). SUEWS has demonstrated good performance against hydrological observations and surface flux observations in several cities in Europe, North America, and Asia (Alexander et al.2016; Ao et al.2018; Havu et al.2022a; Järvi et al.2011; Ward et al.2016). In SUEWS, the seasonal cycle of vegetation phenology is indicated by leaf area index (LAI). Previous studies made in two UK cities and Shanghai, China, have reported that bias in modeled LAI leads to overestimation or underestimation in QE or QH (Ao et al.2018; Ward et al.2016). They highlighted the importance of having an appropriate seasonal cycle of LAI in SUEWS. Omidvar et al. (2022) proposed a workflow to derive LAI-related parameters for SUEWS, but it was intended for fully vegetated areas located mainly on the outskirts of cities. Apart from LAI, the maximum conductance (gmax) is also critical in scaling the surface conductance (gs) and therefore the available energy distribution (Ward et al.2016). However, the impact of LAI-related parameters and gmax on the modeled turbulent fluxes has received insufficient attention in urban areas.

Recently, the module of local-scale CO2 flux (FC) was incorporated into SUEWS (Järvi et al.2019). It was found to give reasonable annual sum, seasonal, and diurnal cycles against observed FC in Helsinki, Finland, suggesting that the bottom-up CO2 emission or uptake estimate in SUEWS can be evaluated by observation-based evidence provided by top-down eddy covariance (EC) measurements. Furthermore, SUEWS shows the potential for broader use, such as quantifying the carbon sequestration potential of urban vegetation (Havu et al.2022a), investigating the spatial variability of CO2 emission, quantifying the contribution of each emission component (Järvi et al.2019), and assisting urban CO2 emission mitigation. However, this module has not yet been evaluated in cities other than Helsinki.

Beijing provides a unique test bed for SUEWS evaluation and application: a megacity with a population of over 21 million and an increasing urbanized area (MHURD2018). The older version of SUEWS (V2017b) has been evaluated and applied in Beijing by Kokkonen et al. (2019), showing good model performance against observed heat fluxes. However, good simulation of turbulent flux does not necessarily imply that the sub-models within give accurate estimates, e.g., LAI and radiative components. Correct presentation of these processes is necessary for more advanced applications like the prediction of surface exchanges of energy under future climate scenarios. Besides, the newly developed FC module has not yet been evaluated in Beijing.

In this paper, we present a comprehensive evaluation of SUEWS V2020b and its ability to simulate surface fluxes of energy and CO2 in Beijing. The main aims of this study are (1) to evaluate the model performance of SUEWS using different (default and site-specific) vegetation parameters against the turbulent flux (QE, QH and FC) measurements and (2) to estimate the contributions of the anthropogenic and biogenic components to the FC with the bottom-up modeling approach of SUEWS. Meanwhile, the impact of the mismatch between the turbulent flux source area and the modeled area is also examined.

2 Model description

SUEWS is an urban land surface model that simulates the surface energy and water balances and CO2 flux at a local (neighborhood) scale (Järvi et al.2011, 2019; Ward et al.2016). In SUEWS, the modeling domain is separated into seven interacting surface types (buildings, paved surfaces, grass, evergreen trees/shrub, deciduous trees/shrubs, bare soil, and waterbody), with a single soil layer below each type. SUEWS is designed to run with surface information (e.g., surface cover fractions) and a minimal amount of model forcing data, including wind speed (U), relative humidity (RH), air temperature (Tair), air pressure (p), precipitation, and incoming solar radiation (Kdown). SUEWS has sub-models for LAI and net all-wave radiation, and users are allowed to modify the parameters of the sub-models based on the information of the modeled domain. In this study, we use SUEWS version V2020b (Havu et al.2022b).

2.1 Leaf area index model

In SUEWS, leaf growth is accumulated when Tair stays above the limit value Tbase,GDD,i consecutively, denoted by growing degree day (GDD). Leaf growth is allowed until GDD reaches the upper boundary GDDfull,i or LAI reaches its maximum (LAImax,i). Similar to the leaf growth period, the leaf-off period is impacted by Tbase,SDD,i, senescence degree day (SDD), and SDDfull,i or LAImin,i. Leaf fall is controlled by Tair or (at high latitudes) by day length (Järvi et al.2014). Here, LAI for vegetation type i at the day of year d (LAId,i) is defined as follows:

(1) LAI d , i = min LAI max , i , LAI d - 1 , i ω 1 , GDD , i GDD d , i ω 2 , GDD , i + LAI d - 1 , i , leaf-on , T base , GDD , i > T d - 1 max LAI min , i , LAI d - 1 , i ω 1 , SDD , i SDD d , i ω 2 , SDD , i + LAI d - 1 , i , leaf-off , T d - 1 < T base , SDD , i ,

where LAImax,i and LAImin,i for each vegetation type can be obtained from literature or determined from observations, ω1/2,GDD/SDD,i represents the growing or senescence rates derived for each study site or uses their default values (Järvi et al.2011; Omidvar et al.2022), and Td−1 is the previous day air temperature mean.

2.2 Radiation fluxes

Kdown is a required variable in the meteorological forcing data, whereas other radiation components are estimated within SUEWS. Outgoing shortwave radiation (Kup) is calculated using a bulk albedo (α) based on the area fraction for each surface type. Incoming longwave radiation (Ldown) is calculated using Tair and RH to estimate the cloud cover and the clear-sky emissivity (Loridan et al.2011), while outgoing longwave radiation (Lup) is estimated by surface emissivity, α, Kdown, Ldown, and Tair (Offerle et al.2003).

2.3 Turbulent heat fluxes

Latent heat flux (QE, W m−2) is calculated using the modified Penman–Monteith equation for each surface type:

(2) Q E = s ( Q N + Q F - Δ Q s ) + ρ c p VPD / r av s + γ ( 1 + r s / r av ) ,

where QN (W m−2) is the net all-wave radiation, QF (W m−2) is the anthropogenic heat flux, ΔQS (W m−2) is the net storage heat flux, ρ (kg m−3) is the air density, cp (Jkg-1K-1) is the specific heat capacity of air at constant pressure, VPD (Pa) is the vapor pressure deficit, s (PaC-1) is the slope of the saturation vapor pressure curve, γ (PaC-1) is the psychrometric constant, rav (s mm−1) is the aerodynamic resistance for water vapor, and rs (s mm−1) is the surface resistance. rs, or its inverse surface conductance gs (mm s−1), has the following form:

(3) g s = 1 r s = i g max , i LAI i LAI max , i fr i × G 1 g ( K down ) g ( Δ q ) g ( T air ) g ( Δ θ ) ,

where gmax,i is the maximum conductance of vegetation type i; fri is the surface fraction of i; G1 is a constant connecting stomatal conductance to canopy conductance; and g(Kdown), gq), g(Tair), and gθ) are environmental response functions on Kdown, specific humidity deficit (Δq), air temperature (Tair), and soil moisture deficit (Δθ), respectively. The functions have the following forms (Ward et al.2016):



(7) T C = ( T H - G 5 ) ( G 5 - T L ) ,


(8) g ( Δ θ ) = 1 - exp ( G 6 ( Δ θ - Δ θ WP ) ) 1 - exp ( - G 6 Δ θ WP ) .

Coefficients G2G6 determine the shape of the curves describing responses of stomatal conductance to each environmental variable. Kdown,max (W m−2) is the maximum incoming solar radiation; TL and TH (C) are the lower and upper limits for temperature at which photosynthesis and transpiration are off, respectively; and ΔθWP (mm) is wilting point deficit. Kdown (W m−2) is model input, Δq (g kg−1) is calculated from model input RH, Tair (C) is either model input or simulated at 2 m height, and Δθ (mm) is simulated within SUEWS (Järvi et al.2017). QH is determined as the residual from the surface energy balance equation:

(9) Q H = Q N + Q F - Δ Q S - Q E .

2.4CO2 flux

The FC module adopts a bottom-up approach to determine the local-scale FC (µmolm-2s-1), accounting for both anthropogenic (FC,ant) and biogenic (FC,bio) components (Järvi et al.2019):

(10) F C = F C , ant + F C , bio = ( F M + F V + F B + F P ) + ( F pho + F res ) .

In Eq. (10), FM is the CO2 emissions from human metabolism, FV is the emissions from traffic, FB is the emissions from buildings (e.g., combustion of natural gas, coal, and wood), FP is the emissions from local-scale point sources, Fpho is the CO2 uptake by photosynthesis, and Fres is the CO2 release by soil and vegetation respiration. Positive values indicate CO2 emissions, and negative values indicate CO2 uptake with respect to the atmosphere. Fpho has a negative sign, while the rest of the FC components have a positive sign.

FC,ant relates to the energy balance through QF. FM and FV are estimated with an inventory approach, i.e., based on population density or traffic rate, and their emission factors (EFs). Hourly CO2 emissions from human metabolism on weekdays or weekends (FM,h,d, µmolm-2s-1) are calculated using

(11) F M , h , d = p d PP h , d AP h , d C M ,

where pd is the daily average population density (cap ha−1), PPh,d is the population diurnal profile by hour, APh,d is the activity level diurnal profile by hour, and CM is the CO2 release per person (µmol CO2 s−1 per capita). The pd and PPh,d reconstruct the diurnal population density cycle. APh,d scales the CM to vary between the nighttime minimum and daytime maximum values (CM(min,max)) to indicate the diurnal cycle of per capita human metabolic intensity.

Hourly traffic CO2 emissions (FV,h,d) on weekdays or weekends are calculated from

(12) F V , h , d = Tr d E c , d H T , d ,

where Trd is the mean daily traffic rate within the study area (vehd-1area-1), HT,d is the diurnal traffic profile, and Ec,d is the traffic EFs for CO2 (kg km−1 per vehicle).

Hourly building CO2 emissions (FB,h,d) on weekdays or weekends are calculated from

(13) F B , h , d = [ fr heat Q F , heat + fr nonheat Q F , base fr QF , base , BEU , d ] × E CO 2  per J ,

where frheat is the fraction of fossil fuels used for heating; QF,heat is the building heat emission at local scale estimated from the heating degree day model (Järvi et al.2011); frnonheat is the fraction of fossil fuels used for energy other than heating (e.g., the use of gas stove for cooking); QF,base is the non-temperature-related anthropogenic heat flux (W m−2) including heat emissions from traffic, human metabolism, and electricity usage; frQF,base,BEU,d is the fraction of the QF,base coming from building energy use on weekdays or weekends; and ECO2 per J is the EF for fuels in building energy use (µmol CO2 J−1). Emissions from single-point sources such as power plants and industrial activities can be included as model inputs.

Fpho relates to the energy balance through LAI and the environmental responses of surface conductance (Eq. 3). Fpho is calculated using

(14) F pho = i ( fr i F pho , max , i LAI i ) × g ( T air ) g ( Δ q ) g ( Δ θ ) g ( K down ) ,

where Fpho,max,i is the maximum photosynthetic rate for vegetation type i.

Soil and vegetation respiration Fres (µmolm-2s-1) follows an exponential dependency on Tair:

(15) F res = i fr i max ( a i exp ( T air b i ) , 0.6 ) ,

where ai and bi are the parameters controlling the temperature dependency and 0.6 µmolm-2s-1 is the minimum respiration in winter time.

Figure 1(a) The location of IAP tower and the land cover type within the 6th Ring Road area of Beijing (MODIS land cover type (MCD12Q1) version 6; Friedl and Sulla-Menashe2019). (b) A satellite image (Google Earth, image ©2022 Maxar Technologies) over the study area. (c) Wind sectors that have been filtered out for data quality control. (d) Urban land use categories (EULUC-China) (Gong et al.2020). (e) Wind direction frequency by season.

3 Study site and measurements

The model domain is a 1 km circle around the 325 m meteorological tower constructed by Institute of Atmospheric Physics, Chinese Academy of Sciences (IAP tower, 3958 N, 11622 E, 60 m above sea level) located in the 6th Ring Road area of Beijing, China (Fig. 1a–d). An EC setup at the height of 47 m on IAP tower continuously measures the surface fluxes of QE, QH and FC using a three-dimensional sonic anemometer (Windmaster, Gill, UK) and an open-path infrared gas analyzer (LI-7500A, LI-COR, USA). In addition, all four radiation components are measured at the height of 140 m using a net radiometer (CNR1, Kipp & Zonen, Netherlands). These measurements are used to evaluate SUEWS model performance. The 1 km radius circle around IAP tower roughly covers 80 % of the accumulated flux footprint area (Liu et al.2012). This area is mainly covered by impervious surfaces (Fig. 1b). Three patches of urban green spaces are situated to the east, south, and west to IAP tower, while the other vegetation is sparsely spread along the roads and near the buildings. Most of the impervious surfaces in the source area are residential buildings (Fig. 1d). A more detailed description of the surroundings and the used instruments can be found in previous publications by Cheng et al. (2018), Liu et al. (2012), and Liu et al. (2021).

The 30 min turbulent flux calculation procedures and the quality controls were described in detail by Cheng et al. (2018). Quality controls such as out-of-limit value removal, spike removal, and a dropout test were conducted on the 10 Hz data during the flux calculation. In order to exclude low-quality data caused by precipitation, dust, or other contamination on the sensor, the records with automatic gain control value  62 were discarded. On top of the procedures by Cheng et al. (2018), the following quality control steps are performed for 30 min turbulent flux observations in the year 2016. (1) Upper- and lower-boundary filtering is performed to remove observations that fall outside specified ranges: QE from 500 to 1000 W m−2, QH from 500 to 1000 W m−2, and FC from 100 to 200 µmolm-2s-1. (2) Spike detection is done to remove flux values outside 3.5 standard deviations of the 3 d moving mean value (Liu et al.2012). (3) Wind direction filtering is performed to remove the wind directions with building heights over 50 m (112–128, 160–243, 314–3) (Kokkonen et al.2019) (Fig. 1c). (4) Stationarity tests are performed to filter out data points with stationary indicator > 30 % (Foken and Wichura1996). The percentage of data removed through these four steps is 0.2 %–0.3 %, 1.7 %–2.7 %, 37.8 %–38.2 %, and 13.1 %–17.4 %, respectively. The numbers of observations retained after quality control are 8017 for QE, 7338 for QH, and 7797 for FC. These 30 min flux observations are resampled to 1 h resolution. In addition, FC is gap-filled using seasonal mean diurnal cycle in order to calculate the seasonal and annual sums (Falge et al.2001).

Wind directions are mainly from the southern and western sectors and the northeastern and eastern sectors before the implementation of the wind direction filtering (Fig. 1 e). In winter, wind from the west is more frequently seen than from the east compared to the other seasons.

To optimize the behavior of LAI, a 6-year time series (2011–2016) of LAI over an adjacent park near the IAP tower is calculated from the atmospherically corrected surface reflectance provided by the USGS Landsat 7 Enhanced Thematic Mapper + (ETM+) (30 m spatial resolution) via the Google Earth Engine Data Catalog (Masek et al.2006). The atmospherically corrected surface reflectance bands have been preprocessed using the scaling factors from the metadata. Next, enhanced vegetation index (EVI) is calculated using the following formula (Huete1997):

(16) EVI = 2.5 × ( NIR - RED ) / ( NIR + 6 × RED - 7.5 × BLUE + 1 ) ,

where NIR, RED, and BLUE are the near-infrared, red, and blue bands, respectively. EVI is further used to calculate LAI with the formula (Boegh et al.2002):

(17) LAI = 3.618 × EVI - 0.118 .

The LAI and air temperature time series are subjected to optimization using the covariance matrix adaptation evolution strategy (CMA-ES) (Appendix A). Before the optimization process, values larger than 10 m2 m−2 and negative values are considered outliers and removed; values during December and January are set to a fixed value, i.e., the average of these months (0.2 m2 m−2), to reduce the noise in winter and improve optimization performance. More details can be found in Appendix A. The related data and code are openly available to reproduce the results (Zheng et al.2022).

4 Model run

4.1 Forcing meteorological data

The reanalysis dataset WFDE5 (Cucchi et al.2021) is used as the forcing data for SUEWS. WFDE5 is a bias-corrected dataset of near-surface meteorological variables specifically suited for land surface modeling. It is derived from the fifth generation of the European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalysis (Hersbach et al.2020). It is provided at 0.5 spatial resolution and at hourly temporal resolution. WFDE5 is evaluated against observed meteorological observations before it is used as the forcing data for SUEWS (Appendix B).

4.2 Land cover

Land cover types and their fractions needed in the model cases are estimated based on aerial images. Paved surfaces account for 46 % of the total area, buildings 24 %, trees/shrub 13 %, grass/lawn 16 %, and water 1 %. The average building height is 19.1 m (Kokkonen et al.2019). According to a field survey conducted in the 6th Ring Road area, the population of deciduous species accounts for 82 % of the total number of woody plants investigated (Ma2019). Therefore, the surface fraction of deciduous trees is set to 11 % and evergreen trees 2 %.

Due to the wind direction filtering (Sect. 3), the actual flux source area “seen” by the EC measurement is biased from the 1 km radius circle around IAP tower. The vegetation fraction is 31 % in the remaining sectors combined, as compared to 29 % in the entire 1 km radius circle (Fig. 1b–c, Table S1 in the Supplement). The model performance in turbulent flux modeling with the land fractions for the remaining sectors is similar to the entire circle (Fig. S1 in the Supplement). Therefore, only the results using the land fractions of the entire circle are demonstrated in the main text.

4.3 Storage heat flux

To calculate the storage heat flux, the Objective Hysteresis Model (OHM) is used (Grimmond and Oke1999). The coefficients for all the surface types follow the previous study by Kokkonen et al. (2019). A large portion of the paved surface is asphalt in the study area. Thus, the coefficients are set to the weighted average values of asphalt surface (AN99) following Ward et al. (2016).

4.4 Human activity

In this study, local parameters of traffic, population dynamics, and building energy use are incorporated in order to estimate QF and CO2 emissions.

Ward et al. (2013)Moriwaki and Kanda (2004)Wen et al. (2020, 2022)Zhang et al. (2014)Wen et al. (2020, 2022)Zhang et al. (2014)Yang et al. (2019)Yang et al. (2019)Cui et al. (2019)MHURD (2018)Zhang et al. (2019)BMBS (2017)Cui et al. (2019)

Table 1Parameters related to anthropogenic heat and CO2 emissions.

Download Print Version | Download XLSX

Figure 2Annual average diurnal cycle of traffic rate for weekday (Trwd) and weekend (Trwe), population density (POP), and the activity profiles for weekdays (APwd) and weekend (APwe) within the 1 km radius circle around IAP tower.


The annual mean weekday and weekend diurnal cycle of traffic rate for each road link in 2017 in the study domain are all extracted from a dataset based on an extensive road traffic monitoring network (Yang et al.2019). For weekends and weekdays, the diurnal traffic cycles are calculated independently. The total hourly traffic rate (vehicle kilometer per hour) is calculated as the sum of the traffic rates, i.e., the product of traffic volume (vehicles per hour) and the road link length (km) from all the road links in the study area. The hourly traffic rates are then summed up to the total daily traffic rates (vehicle kilometer per day) and divided by the total modeling area, yielding Trd. Finally, the diurnal traffic profiles (HT,d) are obtained by normalizing the diurnal cycles of the total hourly traffic rate (Table 1, Fig. 2).

The annual average diurnal cycle of population density within the model domain is obtained from a dataset of hourly population dynamics at 500 m resolution generated from remotely sensed and geospatial data over the years 2015 and 2016 (Zhao et al.2021) (Fig. 2). As there are several residential building areas located within the model domain (Fig. 1d), population density increases in the evening when residents get home from work, remains high throughout the night, and decreases in the morning. Weekdays and weekends share the same diurnal cycle of population density in our study.

Heating in Beijing is dominated by central heating, supplied mainly at the district level. The sources include cogeneration plants fueled by coal or gas and district boilers powered by coal, oil, or gas. Cogeneration plants are usually located at suburban or rural areas, and there are no cogeneration plants within the model domain, meaning that their contribution to CO2 emission is neglected in this study. In comparison, boiler plants are very common: over 5000 coal-fired and 1000 gas-fired heating boilers were located surrounding the populated areas in 2014 (Cui et al.2019). We investigated the boiler plants near the IAP tower through interviews and found that there were at least 11 of them located at multiple directions within a 1.5 km distance from the IAP tower. For the known boiler plants, eight of them have a chimney height lower than 20 m. Thus, their CO2 emissions are very likely to be observed by EC at 47 m during the heating season. Unfortunately, the heating capacity and detailed information regarding the fuel combustion for each boiler plant are either unknown or restricted from access. Therefore, it is challenging to treat the boiler CO2 emissions as point sources. As an alternative, SUEWS first estimates the anthropogenic heat release from heating QF,heat and then converts the heat into local CO2 release using the EF and the fraction of fossil fuels used for heating frheat (Eq. 13).

In 2015, the ratio of district boiler heating capacity to cogeneration plants was 4.2:1 (Zhang et al.2019). Correspondingly, frheat is set to 0.81 to represent the heating capacity from local boiler plants. The value of natural gas consumption over the annual total heating supply is 3.2 × 106tce (tonne coal equivalent) and the value for coal consumption 2.6 × 106tce in 2015. The consumption of the rest of the fuel types is only 3.8 × 105tce (MHURD2018). The EFs of heating supply are 96.51 Tg CO2 (1019 J)−1 for coal-fired boilers,and 56.17 Tg CO2 (1019 J)−1 for gas-fired boilers (Du et al.2018). Therefore, the EF for fuels used in building energy use (ECO2 per J) in SUEWS takes the average of the EFs of natural gas and coal weighted by their annual consumption, i.e., 0.1688 µmol CO2 J−1 (Table 1). In addition, SUEWS needs a temperature limit (base temperature, Tbase_HC) to indicate when heating takes place in the heating degree day model. Central heating usually starts around 15 November and lasts until 15 March, when the outdoor air temperature is around 12 C. Therefore, this value is given to SUEWS Tbase_HC to replace the default value (18.2 C) (Järvi et al.2011).

Statistics showed that urban household living consumed liquefied petroleum gas at a rate of 27.9 × 107kgce (kilogram of coal equivalent), gas 17.1 × 108kgce, and electricity 20.46 × 108kgce in 2016 (BMBS2017), indicating that 50 % of the household energy use involved on-site CO2 emissions. Therefore, the non-heating fraction (frnonheat) is set to 0.5.

4.5 Evaluation design

Two different groups of SUEWS runs are made around the IAP tower (Fig. 1c). The first run over the years 2009 to 2011 is to evaluate the modeled radiation components against observations. The first 16 months are the spin-up period, and the actual model evaluation is made from May 2010 to June 2011, when the radiation observations are available. The second SUEWS run from 2015–2016 is to evaluate the turbulent fluxes. The first year is used as a spin-up period, and only the year 2016 is used in the evaluation.

The model performance of radiation fluxes is evaluated prior to the simulation of turbulent heat fluxes. The results show that SUEWS is applicable to provide realistic estimates of radiation fluxes in the study area despite the absence of site-specific parameters (Appendix C). The calculation of radiation fluxes is mostly dependent on the land cover fractions under the current scheme adopted by SUEWS (Loridan et al.2011; Offerle et al.2003). No visible change in land use type is observed according to satellite images from Google Earth within the modeled area between 2010 and 2016 (figure not shown). Therefore, we assume that the evaluation of radiation fluxes using observations in the years 2011–2012 holds true in 2016.

Järvi et al. (2011)Järvi et al. (2011)Havu et al. (2022a)Havu et al. (2022a)Havu et al. (2022a)Havu et al. (2022a)Havu et al. (2022a)Havu et al. (2022a)Havu et al. (2022a)Havu et al. (2022a)Järvi et al. (2014)Ward et al. (2016)Ward et al. (2016)Havu et al. (2022a)Järvi et al. (2019)Havu et al. (2022a)Järvi et al. (2019)Havu et al. (2022a)Kokkonen et al. (2019)

Table 2SUEWS biogenic model parameters used for all case runs in this study.

a Evergreen tree. b Deciduous tree.

Download Print Version | Download XLSX

4.5.1 Sensitivity to vegetation-related parameters

In order to test the model's sensitivity of radiation and turbulent fluxes to vegetation-related parameters, four model cases are designed as follows.

  1. The “Base” case. Control run where model parameters are considered “default” following Kokkonen et al. (2019) (Tables S2–S4). The exceptions are (1) parameters in the environmental response functions (Eqs. 47) of surface conductance that follow those by Havu et al. (2022a) (Table 2) because the product of response functions calculated following Kokkonen et al. (2019) is too low (95th percentile = 0.19) to obtain a realistic estimate of photosynthetic rate and (2) biogenic parameters for Eqs. (14) and (15) and soil properties that were updated following Havu et al. (2022a) (Table 2). In addition, Fpho,max for grass and lawn (5.5 µmolm-2s-1) is obtained from EC observations of CO2 flux made over urban lawn in Helsinki in summer 2021 by fitting the conductance parameters (Fpho,max,grass, G2G6) to the observations following Järvi et al. (2019) (Appendix D).

  2. The “LAI” case. The same as the base case but the parameters for Eq. (1) describing the annual behavior of LAI are optimized using remotely sensed LAI and CMA-ES (Appendix A). The new optimized LAI parameters are compared with the base case in Table 3.

  3. The “gs” case. The same as the base case but with gmax values collected from observational studies over vegetation species in Beijing (Appendix E). The site-specific gmax values are in general lower than the values used by the base case (Table 4).

  4. The “gs_LAI” case. The same as the base case but with both LAI and gmax modified as described in the LAI and gs cases.

Table 3Comparison in leaf area index (LAI) parameters between the base case and LAI and gs_LAI cases. All vegetation types (evergreen tree, deciduous tree, and grass) use the same LAI parameters within one case. The base case values are the same as in Järvi et al. (2011).

Download Print Version | Download XLSX

4.5.2 Sensitivity to radius of modeled area

SUEWS output will be evaluated against EC-measured turbulent fluxes, but the challenge is that the source area of the observations is different to the exact modeling domain. To consider the impact of the chosen modeling domain on model evaluation, we designed three additional model cases where circular areas around the IAP tower with different radii are considered. The default run is with the 1 km radius circle, but SUEWS is also run within 500, 750, and 1500 m circular areas, corresponding to the flux footprint of roughly 60 %, 70 %, and 80 %–90 %, respectively (Liu et al.2012). Used vegetation parameters are as in the gs_LAI case, but land surface fractions, population, and traffic parameters are modified accordingly (Appendix F).

4.6 Statistical metrics for model evaluation

Common statistical metrics are used to quantify the model performance, including the coefficient of determination (R2), root-mean-square error (RMSE), and mean bias error (MBE). Simple linear regression is used to estimate the relationship between the model output and the observations, and the square of the correlation coefficient is taken as R2. The other statistical metrics are defined as follows:


where yi^ is the modeled value and yi the measured value. Statistical metrics are calculated at annual and seasonal scale, i.e., DJF (winter), MAM (spring), JJA (summer), and SON (autumn).

Table 4Comparison of maximum conductances (gmax) between the base case and gs and gs_LAI cases. Parameters for the base case follow Järvi et al. (2011). More details can be found in Appendix E.

Download Print Version | Download XLSX

5 Results and discussion

5.1 Seasonal dynamics of the optimized LAI

The control base case simulates the onset of leaf growth and the ending of senescence reasonably well (Fig. 3). The performance of LAI modeling is further improved after the optimization (Appendix A). In the base case, the modeled LAI starts to increase rapidly from day-of-year (DOY) 70 and plateaus at DOY 105, which is too early when compared to the remotely sensed LAI (Landsat 7 LAI). The optimized LAI starts to grow at the same time but slightly slower and peaks 20 d later than in the base case. In autumn, LAI modeled by the base case drops rapidly at DOY 310, while the optimized LAI starts to decline rapidly at DOY 267. The LAI model with optimized parameters is better at capturing the behavior of senescence than in the base case.

Figure 3Normalized LAI in 2016, where the value 0 (1) represents the minimum (maximum) of LAI. “Base” denotes the modeled LAI for the base case, “Landsat” denotes the remotely sensed LAI of the adjacent park near the IAP tower, and “optimized” denotes the modeled LAI using the parameters derived with CMA-ES (Appendix A).


Although previous studies suggested that LAI was generally modeled well using default parameters following Järvi et al. (2011), Ward et al. (2016) reported that leaf growth is reached too early and suddenly in spring in two UK cities. In contrast, Ao et al. (2018) showed that the simulated LAI might be lower than reality when vegetation remained green in winter and spring in Shanghai, China. These lead to the bias in gs and therefore QE estimates (Eqs. 2 and 3). In Beijing, the rainy season lasts from May to October, while the rest of the year is the dry season (Liu et al.2012). It is possible that the distinct dry season leads to a lack of soil moisture in spring and autumn and thus influences the LAI seasonal dynamics if there is no external water input (Omidvar et al.2022). However, the urban green spaces in Beijing are usually sufficiently or even excessively irrigated (Zhang et al.2017). Observations also provided evidence to support the relationship between air temperature and phenological dynamics in the urban environment in Beijing (Lu et al.2006; Luo et al.2007). Therefore, the air temperature-dependent LAI model is applicable in Beijing, but the “default” LAI parameters might be not suitable. We recommend evaluating the LAI model when SUEWS is applied to a different city and deriving the optimal LAI parameters if necessary.

Figure 4Annual and seasonal mean diurnal cycles of observed and modeled (a–e) latent heat flux (QE) and (f–j) sensible heat flux (QH) for the four model cases (base, gs, LAI, and gs_LAI) in the year 2016. The shaded area denotes the interquartile range.


Figure 5Model performance statistics R2, RMSE, and MBE of (a–c) latent heat flux (QE) and (d–f) sensible heat flux (QH) for the four model cases (base, gs, LAI, and gs_LAI) in the year 2016.


5.2 Evaluation of turbulent heat flux modeling

Both observed QE and QH reach their maxima around noon (Fig. 4). The observed QE has the largest amplitude during the summer months, while QH during the spring months. All four model cases capture their diurnal cycles, but large differences in the amplitude and model performance are observed among the model cases (Fig. 5).

In the base case, the model overestimates QE (with MBE from 7.4 to 48.6 W m−2) except in winter months (Fig. 5). With the optimized LAI (the LAI case), model performance in QE remains virtually unchanged, with RMSE (12.1–94.1 W m−2) and R2 (0.17–0.53) when compared to the base case (with RMSE 11.7–96.1 W m−2 and R2 0.20–0.51) (Fig. 5a–c). QE is to a large extent determined by surface conductance, which is scaled by gmax of each vegetated surface for the modeled area (Eq. 3). Clear improvement is observed when the local gmax is used (the gs case), especially during summer. The overestimation of QE is largely reduced, RMSE drops to 11.3–54.7 W m−2, and R2 increases to 0.25–0.61. With both the optimized LAI and the local gmax introduced (the gs_LAI case), the R2 is similar to the gs case, while RMSE slightly decreases by 1–2 W m−2 in spring and summer compared to the gs case. QE is underestimated throughout the day in winter (MBE=8.3 W m−2). Observational studies have shown that combustion-derived water vapor often contributes 5 %–10 % of total urban humidity during the heating season (Fiorella et al.2018; Liu et al.2022; Salmon et al.2017). The observed QE ranges from 0 to 30 W m−2 in winter when vegetation is dormant; this suggests that combustion and evaporation, as the dominant sources of QE, might lead to QE at this magnitude. The anthropogenic water vapor release might be underestimated by SUEWS in winter.

In SUEWS, QH is estimated as the residual of energy balance and is therefore directly affected by the modeled QE. As a result of overestimating QE in the base case, QH is greatly underestimated with MBE of 51.1–14.4 W m−2 and RMSE of 60.9–116.0 W m−2. R2 values in summer months and autumn months are lower than 0.1. The model performance is barely improved by the optimized LAI (the LAI case) but it is noticeably improved after the local gmax is introduced (the gs case) (Fig. 5d–f). The best model performance for QH is obtained by the gs_LAI case, decreasing the magnitudes of MBE to 13.6–15.2 W m−2 and RMSE to 60.1–75.3 W m−2 and increasing R2 to 0.40–0.64 (Fig. 5d–e).

QH is also influenced by QF. Nighttime QF in summer might be overestimated, leading to the overestimation in QH. Turbulent heat fluxes are also related to ΔQS. Both QE and QH correlate negatively with ΔQS (Eqs. 2 and 9). For instance, QE is underestimated, while QH is overestimated at noon in summer in the gs_LAI case. The decrease in ΔQS will lead to a simultaneous increase in QE and QH, lowering the bias of QH while increasing the bias of QE. Therefore, the adjustment of ΔQS can hardly improve the QE and QH modeling at the same time.

Our results suggest that the model performance of heat fluxes is more sensitive to the adjustment of gmax than to the change in LAI seasonal dynamics. By incorporating local LAI and gmax, SUEWS simulates the heat fluxes noticeably better, increasing R2 by 0.03 (0.30) and decreasing RMSE by 27.0 (23.7) W m−2 for QE (QH) compared to the base case.

Figure 6Annual and seasonal average diurnal cycles of observed and modeled CO2 flux (FC) for the four model cases (base, gs, LAI, and gs_LAI) in the year 2016. The shaded area denotes the interquartile range.


Figure 7Seasonal average diurnal cycles of modeled CO2 flux (FC) components in the gs_LAI case in 2016. FV denotes the CO2 emissions from on-road traffic, FB denotes buildings, FM denotes human metabolism, Fres denotes vegetation and soil respiration, and Fpho denotes the CO2 uptake by vegetation photosynthesis. Positive values indicate emissions of CO2, and negative values indicate the uptake of CO2 with respect to the atmosphere. The shaded area denotes the interquartile range.


5.3 Evaluation of CO2 flux modeling

5.3.1 Model performance

SUEWS basically reproduces the average annual and seasonal diurnal cycle of observed FC (Fig. 6). The diurnal behavior is dominated by on-road traffic emission, which reaches 22.6 and 23.0 µmolm-2s-1 for the morning peak and afternoon peak, respectively, during the rush hours (Fig. 7). Human metabolism (maximum 4.8 µmolm-2s-1) is the second largest source of CO2 emissions. In winter, the building CO2 emission has a maximum of 6.7 µmolm-2s-1 in the daytime. The maximum photosynthesis rate is 5.9 µmolm-2s-1 around noon in summer, while soil and vegetation respiration constantly serves as a CO2 source with a rate lower than 3.6 µmolm-2s-1.

Figure 8Model performance statistics for (a) R2, (b) RMSE, and (c) MBE of CO2 flux (FC) for the four cases (base, gs, LAI, and gs_LAI) in the year 2016.


Each of the model performance statistics of FC is of a similar magnitude among cases, indicating the FC modeling is far less sensitive to the choices of gmax and LAI-related parameters than QE and QH shown in Sect. 5.2 (Fig. 8). Under the current parameterizations, Fres considers only air temperature (Eq. 15). The adjustments of gmax and LAI parameters affect the modeled heat fluxes, influencing 2 m air temperature, and finally Fres, but the difference in annual CO2 release from respiration is less than 0.01 kgCm-2yr-1 among cases. Fpho is sensitive to the adjustments of gmax and LAI parameters. In the base case, the large values of gmax allow relatively high evapotranspiration (namely QE). As a result, the average Δθ during January and June is larger than 105 mm, which is only 27 mm lower than the wilting point deficit (ΔθWP). The dry soil lowers the surface conductance and photosynthetic CO2 uptake through the limiting function of gθ) (Eq. 8). As the local gmax is introduced, the soil remains moister with Δθ lower than 75 mm throughout the year, allowing more favorable conditions for the photosynthetic CO2 assimilation. The CO2 assimilated through photosynthesis is 0.57 kgCm-2yr-1 in the gs case, which is 0.21 kgCm-2yr-1 higher than in the base case. The LAI reduction in spring and autumn in the gs_LAI case, on the other hand, directly limits surface conductance and photosynthesis (Eq. 14), leading to a decrease of 0.07 kgCm-2yr-1 in annual photosynthetic CO2 uptake when compared to the gs case.

In SUEWS, photosynthetic and respiration rates are proportional to the fractions of vegetated surfaces, which account for only 29 % of the modeled area. The magnitude of Fpho is substantially lower than the traffic emission, making the effect of photosynthesis, as well as its response to the adjustments of gmax and LAI parameters, hardly visible in the FC diurnal cycles (Fig. 6).

Expectedly, SUEWS has difficulty in capturing the hourly variability of FC, resulting in the overall low R2 (0.16–0.22) and high RMSE (12.9–16.4 µmolm-2s-1) (Fig. 8a). On the one hand, observed FC has great variability at the hourly scale throughout the year, as indicated by the large interquartile range (Fig. 6). On the other hand, under the current parameterization, two of the anthropogenic FC components are static: modeled traffic emission diurnal cycle is only dependent on whether it is a weekend or a weekday, and modeled human metabolism diurnal cycle is invariable throughout the year (Fig. 7), making it difficult to capture the extreme values. Other urban FC bottom-up modeling studies also reported similar challenges in modeling FC hourly variability in Helsinki and in Basel (Järvi et al.2019; Stagakis et al.2022).

SUEWS gives a reasonable estimate of annual accumulated FC (8.6 kgCm-2yr-1), which is 15 % higher than the observed gap-filled value (7.5 kgCm-2yr-1). FC is overestimated in all seasons, with the lowest MBE (0.2 µmolm-2s-1) in winter and the highest MBE (5.7 µmolm-2s-1) in summer (Fig. 8c).

There are multiple reasons to explain the difficulty in accurately capturing the diurnal cycle of the observed FC for each season. First, the underlying seasonal variation in the diurnal wind pattern makes the FC diurnal cycle from the NW quadrant more “seen” in winter and spring, while the diurnal cycle from the SE quadrant is more seen in summer and autumn (Fig. S2a–d). Second, the observed FC varies noticeably with wind direction (Fig. S2e–h). An evident afternoon and evening FC peak is observed from the NE and SE throughout the year. This might be attributed to the relatively high traffic volumes and more severe traffic congestion in the afternoon and evening than in the morning, especially at the nearby crossroads (Fig. 1d). The morning FC peak only comes from the NE and NW and only during winter months, suggesting that there might be seasonal variation in the traffic pattern that is not captured by SUEWS. The noticeable increase in FC in winter from the SE and NW indicates that the building heating emissions might mainly originate from these two sectors. However, SUEWS cannot consider this FC's wind direction dependency as it simulates the overall flux from the simulation domain.

Apart from the wind direction, atmospheric stability influences the real-time footprint fetch of FC (Crawford and Christen2015), but this is not considered in our study. Furthermore, there might be biases in simulating the seasonal cycles of FC components. SUEWS might underestimate the vegetation photosynthetic rate or overestimate the CO2 release from respiration due to the lack of site-specific biogenic parameters. Nonetheless, the model performance over the FC diurnal cycle is reasonably good compared to a previous study (Järvi et al.2019).

Uncertainties in traffic emission originate from the traffic rates and EFs. SUEWS adopts static traffic EFs and neglects the relationship between traffic emission and Tair as reported by Alvarez and Weilenmann (2012) and Fontaras et al. (2017). In order to examine the impact of seasonal variation of Tair in traffic emission, correction is conducted using the regression function following Zhang et al. (2021), but only a marginal difference is seen at monthly scale: a difference of 3 % in winter, 2 % in spring,  0 % in summer, and 1 % in autumn. Therefore, we believe that the static traffic EFs adopted by SUEWS can provide reasonable traffic emission values without considering the seasonal dynamics of Tair.

Järvi et al. (2019) reported that using a different coefficient of CO2 release per capita (CM) led to a 6 % decrease in the human metabolic CO2 emission estimate. If CM is set to a daily mean value of 242 µmolm-2s-1 (Prairie and Duarte2007) instead of the current values (Table 1), the human metabolic emission will increase and the annual FC will be 4 % higher than the original estimate.

Building emissions are calculated based on the QF estimates and heating fraction. Modeled average QF in December is 52.7 W m−2, which is higher than another model estimate (21.6 W m−2) in the modeled area (Wang et al.2020). Observations of QF are rarely available, and thus these QF estimates have not yet been validated. The representativeness of the heating fraction estimated from yearbook statistics is yet to be examined because the location and heating capacity of heating boilers within the modeled area is unknown. However, the building emission estimate (0.97 kgCm-2yr-1) falls in the range of estimates ( 0–3.0 kgCm-2yr-1) by other cities (Björkegren and Grimmond2018; Christen et al.2011; Järvi et al.2019; Moriwaki and Kanda2004).

The modeled CO2 release by respiration is larger than CO2 assimilated through photosynthesis in our study. At the annual scale, urban vegetative surfaces can have a net CO2 uptake (Awal et al.2010; Konopka et al.2021) but may also have a net emission (Peters and McFadden2012). Admittedly, bias might exist in biogenic CO2 flux estimates since the parameters used in this study are derived from the observations over street trees in Helsinki and over a lawn at Ossinlampi, Finland, where the climate and vegetative species are different from Beijing. With these parameters, the model might underestimate the CO2 sequestrated by the local vegetation or overestimate the CO2 release by respiration.

Figure 9(a) Land use fraction and mean daily traffic rate (Trwd/we) and (b) accumulated CO2 flux (FC) for the modeling domain, with the radius ranging from 500 m to 1500 m. “Wd” denotes weekday, “we” denotes weekend, “Mod” denotes modeled FC, and “Obs” denotes observed FC. Note that in (b) the lines for Mod 500 m and Mod 750 m nearly overlap, while the lines for Mod 1000 m and Mod 1500 m nearly overlap.


5.3.2 The impact of the modeling domain size

The surroundings of the IAP tower are heterogeneous in terms of land surface fraction and mean daily traffic rate (Fig. 9a). The fraction of vegetated surfaces is higher closer to the tower than further away due to the green spaces adjoining the IAP tower (Fig. 1b). Additionally, there is a traffic hot spot on the 3rd Ring North Road located 850 m to the south of IAP tower (Fig. 1d), where the traffic rate is 2 to 7 times the value for the other roads inside the circle of a 1000 m radius (figure not shown). A large increase of 26 % in daily traffic rate is seen when the radius of modeling domain is 1000 or 1500 m when compared to domains with lower radii (Fig. 9a). Thus, the modeled annual accumulated FC largely depends on the modeling domain size chosen, giving estimates of 7.4, 7.6, 8.6 ,and 8.7 kgCm-2yr-1 for the radii of 500, 750, 1000, and 1500 m, respectively. Observational annual FC (7.5 kgCm-2yr-1) falls within this range, which indicates the good model performance of SUEWS (Fig. 9b).

The turbulent flux modeling is usually evaluated over a fixed extent, such as a circle with a certain radius, to approximate the flux source area (Demuzere et al.2017; Järvi et al.2019). However, when a circle with the radius 1000 m is selected to approximate the  80 % footprint fetch in our study, SUEWS does not give the closest estimate of annual FC. This can be explained by the mismatch between the modeling domain and the real flux source area – the single fixed-extent modeled area cannot perfectly represent the land surface characteristics (e.g., the nonuniform land cover and human activities), biassing turbulent flux modeling (Chu et al.2021; Laine et al.2009). First, the accumulated footprint area of the observed fluxes is irregular in shape and varies with time (Liu et al.2012). Second, the relative contribution to flux from the land surface decreases as the distance to the measurement instrument increases (Christen et al.2011; Rebmann et al.2005). Thus, when the modeling domain is a 1000 m radius circle, the model might underestimate the relative contribution from the adjacent vegetated surface and overestimate the contribution of the traffic hot spot at the edge of 80 % footprint fetch.

Regardless of the modeling domain size, traffic is the dominant CO2 source, contributing 59 %–70 % to the total CO2 emissions, followed by human metabolism (14 %–18 %), building energy use (11 %–14 %), and CO2 release by vegetation and soil respiration (6 %–10 %). Vegetation photosynthesis offsets only 5 %–10 % of the total annual CO2 emissions (Fig. 10). Several bottom-up modeling studies show that the on-road traffic is the greatest source in a densely built neighborhood, contributing to 70 % in central London, 61 % in Helsinki, 53 %–78 % in Tokyo, and 70 % in Vancouver, while human metabolism also plays an important role, contributing 5 %–39 % to the annual total FC (Björkegren and Grimmond2018; Christen et al.2011; Järvi et al.2019; Moriwaki and Kanda2004). Our results are in general agreement with these studies. The contribution of the local building emissions within the study area is more variable among cities: a contribution of 70 % is reported in Basel (Stagakis et al.2022), while  0 % is reported in Helsinki (Järvi et al.2019), and our estimate falls in this range. The direct CO2 sequestration by urban plants is minor compared with the total CO2 emissions in this densely built neighborhood, which is in general agreement with Pataki et al. (2011) and Christen et al. (2011). For more accurate biogenic component estimates in Beijing, photosynthetic and respiration observations over local species are needed in the future.

Figure 10The contribution to total CO2 emissions by each component at the annual scale for the modeled area with a radius ranging from 500 to 1500 m. Note that for photosynthesis the percentages denote the offset against total CO2 emissions.


6 Conclusions

A correct description of vegetation is vital in order to simulate the energy and CO2 fluxes over urban surfaces using the urban land surface model SUEWS. In this study, the impact of selecting appropriate vegetation parameters, including LAI parameters and gmax, on simulating the surface fluxes is examined in Beijing, China. In addition, the newly developed CO2 emissions module in SUEWS is evaluated against EC measurements.

The model performance of heat fluxes (QE and QH) is more sensitive to the adjustment of gmax than to the change in LAI seasonal dynamics in our study area. The LAI model has been improved by using CMA-ES to optimize the LAI parameters with a remotely sensed LAI product, providing more realistic vegetation phenological dynamics, especially for the senescence season. The parameter gmax was parameterized according to leaf-level observations of vegetation species in Beijing. By incorporating local LAI and gmax, SUEWS simulates the heat fluxes noticeably better, increasing R2 by 0.03 (0.30) and decreasing RMSE by 27.0 (23.7) W m−2 for QE (QH), and showing more realistic seasonal dynamics when compared to EC observations.

FC modeling appears to be less sensitive to the choice of LAI-related parameters and gmax. Only one of the FC components, Fpho, responds noticeably to them. SUEWS can catch the general diurnal and seasonal behavior of FC but tends to overestimate FC, especially in summer months. We also tested the influence of chosen modeling domain size on simulated FC. By selecting the modeled radii of circular areas ranging from 500 to 1500 m (i.e., an accumulated footprint area from 60 % to 80 %–90 %), the modeled annual FC ranges from 7.4 to 8.7 kgCm-2yr-1, which is comparable with the EC observations (7.5 kgCm-2yr-1). This shows that the model also performs well on the annual scale. Regardless of the modeling domain size, traffic is the dominant CO2 source, contributing 59 %–70 % to the total CO2 emissions, followed by human metabolism (14 %–18 %), buildings (11 %–14 %), and CO2 release by vegetation and soil respiration (6 %–10 %). Vegetation photosynthesis offsets only 5 %–10 % of the CO2 emissions.

We highlight the importance of choosing more site-specific LAI parameters or gmax when SUEWS is used for heat flux modeling before the more advanced application such as urban climate and hydrological modeling. Observations are needed to support more accurate parameterizations of biogenic CO2 fluxes. We believe that the bottom-up approach to model FC by SUEWS can be a promising tool in capturing the CO2 emission hot spots, quantifying the relative contribution of the local CO2 sources, and assisting us in mitigating urban CO2 emissions.

Appendix A: Workflow of LAI parameters optimization

In SUEWS, LAI influences the surface conductance and subsequently QE and Fpho (Sect. 2). A workflow for parameter derivation for the LAI sub-model based on remotely sensed data is designed for natural ecosystems (Omidvar et al.2022). However, vegetation in urban areas behaves differently from natural ecosystems (Zhang et al.2022) and needs to be considered separately. Therefore, we propose a workflow to obtain the parameters for urban area based on remotely sensed LAI and covariance matrix adaptation evolution strategy (CMA-ES). This workflow can also be applied to natural ecosystems. The related data and codes are openly available (Zheng et al.2022).

Covariance matric adaptation evolution strategy (CMA-ES) is one of the strategies for numerical optimization of non-convex problems. It is based on the principle of biological evolution. The evolution strategy takes a certain number of individuals (candidate solutions) in a stochastic way, selects individuals based on the fitness, and repeats this process for generations so that a better or an optimal solution is obtained. Adaptation of the covariance matrix amounts to learning a second-order model of the underlying objective function. Compared with classic optimization methods, CMA-ES requires neither derivatives nor an objective function; it only requires the ranking of candidate solutions. Besides, CMA-ES outranks many of other optimization algorithms, performing especially strongly for “difficult functions” or larger dimensional search spaces (Hansen et al.2010).

Figure A1(a) Time series of the input LAI, the interpolated LAI, and the predicted LAI and a (b) comparison of the predicted LAI against the observed LAI.


Taking our study area as an example, the LAI parameters are optimized as follows:

  1. LAI derivation. A 6-year time series (2011–2016) of LAI of the vegetation in a park near the IAP tower is calculated from the atmospherically corrected surface reflectance provided by USGS Landsat 7 Enhanced Thematic Mapper + (ETM+) (30 m spatial resolution) via the Google Earth Engine Data Catalog (Masek et al.2006). The time series is treated as the “original LAI”.

  2. Spikes removal. There are outliers in the LAI time series caused by instrument problems, uncertainties in the retrieval algorithm, and cloud contamination. Values larger than 10 m2 m−2 and negative values are first removed. The LAI values during December and January are set to a fixed value, i.e., the average of these months (0.2 m2 m−2), in order to reduce the noise in winter and improve the optimization performance.

  3. Scaling the original LAI to the canopy level. The original LAI might be noticeably lower than the measured LAI at the canopy level over a homogeneous vegetated surface. Nonetheless, the original LAI provides the signals of vegetation phenology (e.g., leaf out, peak growing season, leaf fall). In order to give a more realistic estimate of LAI at the canopy level, the original LAI needs to be scaled. Therefore, LAImax and LAImin need to be given manually, preferably based on observational studies over local species. Here, the original LAI is scaled to allow the optimized LAI to reach 5–6 m2 m−2 in the peak growing season as reported by an observational study in Beijing (Wang et al.2021). The canopy-level LAI is marked as the “input LAI” for the process of optimization and marked as the “observed LAI” for the process of evaluation.

  4. Interpolation. The observed LAI values are linearly interpolated between values to obtain a daily time series, marked as the “interpolated LAI”.

  5. Parameter derivation using CMA-ES. The time series of interpolated LAI and Tair are subjected to CMA-ES to optimize the parameters. Using the LAI model and the derived parameters, LAI is calculated and marked as the “predicted LAI”.

Despite the spike removal step, noise is observed in the LAI time series such as the abrupt drop and recovery during summer. Smoothing techniques are used to extract valid information from the satellite data with noise. In other words, the overall pattern (e.g., the moving average) of the time series is the information of interest, and CMA-ES has successfully reproduced the overall pattern out of the input LAI time series contaminated by noise (Fig. A1a). Moreover, the LAI has been “scaled” as shown in the scaling step to allow the output LAI to reach 5–6 m2 m−2 in order to avoid the underestimation of the predicted values. Therefore, the noise in the input LAI during summer is kept since it has little impact on the outcome. We conclude that the model performance is overall good (with R2= 0.74 and RMSE= 1.2 m2 m−2) (Fig. A1b).

Appendix B: Evaluation of WFDE5 reanalysis against observed meteorological variables

To force SUEWS, local-scale meteorological data within the inertial sublayer is required. However, they can be unavailable for the area and period desired. Reanalyses provide spatially and temporally complete datasets, which might make the modeling run easier for users. Kokkonen et al. (2018) and Kokkonen et al. (2019) evaluated one of the reanalyses, WATCH Forcing Data ERA-Interim (WFDEI), suggesting that WFDEI can serve as the forcing of SUEWS but should be corrected beforehand when the bias is large.

WFDE5 is a bias-corrected dataset of near-surface meteorological variables derived from the fifth generation of the European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalysis (ERA5) (Cucchi et al.2021). It is generated using the same methodology as WFDEI, and provides a single layer at 0.5 spatial resolution and hourly temporal resolution. The evaluation of WFDE5 and the use of it as forcing data to SUEWS have been neglected so far. Here, we compare WFDE5 against observed meteorological variables including air temperature (Tair), relative humidity (RH), wind velocity (U) at 47 m, and incoming shortwave radiation (Kdown) at 140 m on the IAP tower (Liu et al.2012). The evaluation of Kdown is conducted from May 2010 to June 2011, and the rest is conducted from January 2010 to December 2011. All the observed variables are resampled from 30 min to 1 h resolution.

With the difference in height for each meteorological variable (Table B1), however, WFDE5 is close to the observed data as a whole (Fig. B1). Compared with the observed data, WFDE5 Tair is lower, RH is higher, U is lower, and Kdown is higher. WFDE5 may underestimate Tair and overestimate RH due to it neglecting the urban anthropogenic heat release; WFDE5 might overestimate Kdown due to insufficient consideration of aerosol's effect in decreasing solar radiation received by the urban surface. The lower U of WFDE5 can be explained by the lower height compared with the observations (Table B1). If the WFDE5 U is adopted as a forcing of SUEWS, aerodynamic resistance might be overestimated, and therefore QE might be underestimated.

Table B1The height of WFDE5 and observed meteorological variables.

Download Print Version | Download XLSX

Admittedly, some meteorological variables of WFDE5 correlate poorly with the observations in a particular season (e.g., R2= 0.13 for U in JJA). However, the overall high R2 and low RMSE and MBE in magnitude suggest that WFDE5 provides reasonably good estimates of each meteorological variable (Table B2). Therefore, WFDE5 is adopted as the forcing data of SUEWS in this study.

Figure B1Comparison of WFDE5 and observed meteorological variables including (a) air temperature (Tair), (b) relative humidity (RH), (c) wind velocity (U), and (d) incoming solar radiation (Kdown) at an hourly resolution.


Table B2Statistics for WFDE5 compared against the observed meteorological variables.

Download Print Version | Download XLSX

Appendix C: Evaluation of radiation fluxes

The radiation parameterization scheme Net All-wave Radiation Parameterization (NARP) does not involve gmax or LAI. As expected, the four experiments (Sect. 4.5.1) give identical radiation flux components in the output. Therefore, only the gs_LAI case is further analyzed here.

The model performance of SUEWS in simulating radiation fluxes is good (Fig. C1, Table C1), and it can reproduce the diurnal cycle of each radiation flux well (Fig. C2). Kup is overestimated in all seasons, with MBE ranging from 0 to 10 W m−2. This overestimation can at least partly be explained by the positive bias (MBE > 15 W m−2 in all seasons) of WFDE5 Kdown when compared to the observed Kdown (Appendix B). The overestimation in Kup can also be caused by surface albedo. Observational studies have shown that the urban surface albedo near the IAP tower varies between 0.1 and 0.15 with season (Jiang et al.2007; Miao et al.2012). The annual bulk albedo for the modeling domain given to SUEWS is 0.14, which is relatively high but still consistent with the observations. A larger positive bias in Kup is observed in summer than in winter. Surface albedo is influenced by many factors, such as surface wetness and street canyon trapping effect (Ao et al.2016; Dou et al.2019), which have not yet been considered by SUEWS. By simply (1) adjusting the albedos for surface types following Ward et al. (2016) and (2) allowing albedo for vegetation to vary from a lower value in summer to a higher value in winter (Table S5), the RMSE for Kup decreases for all seasons, especially in summer (from 18.0 to 13.6 W m−2), but this only has a minor impact on QN modeling (Table S6).

Table C1SUEWS model performance statistics for radiation fluxes. The abbreviations are the same as those in Fig. C1. Note that Kdown is an input of SUEWS, while the others are outputs of SUEWS.

Download Print Version | Download XLSX

The average seasonal and diurnal cycles of Ldown are well captured by the model (Fig. C2i–l), although R2 (0.68–0.88) is lower than with other radiation fluxes. The difference might result from the discrepancy between the observed and modeled cloud fraction as clear skies and overcast conditions are especially difficult to capture using the radiation model NARP, as also reported by Ao et al. (2016) and Ward et al. (2016). The diurnal amplitude of Lup is slightly overestimated by SUEWS (Fig. C2m–p). The model tends to overestimate Lup, particularly at high values of Lup (Fig. C1m–p). The values of emissivity of building materials used in SUEWS might be slightly higher than in reality. Lup is also dependent on Kdown in NARP. Therefore, the overestimation of Lup can partly be explained by the overestimated Kdown provided by WFDE5, especially around noon and in summer.

We conclude that SUEWS is applicable to provide a realistic estimate of radiation fluxes in Beijing, in general accordance with previous studies (Järvi et al.2014; Karsisto et al.2016; Ward et al.2016), despite the absence of site-specific parameters.

Figure C1Input or modeled vs. observed hourly radiation fluxes, including (a–d) incoming solar radiation (Kdown), (e–h) outgoing shortwave radiation (Kup), (i–l) incoming longwave radiation (Lup), (m–p) outgoing longwave radiation (Lup), and (q–t) net radiation (QN) from May 2010 to June 2011. Note that while Kdown is input, the rest are model output.


Figure C2Average diurnal cycle of input or modeled and observed hourly radiation fluxes by season, including (a–d) incoming solar radiation (Kdown), (e–h) outgoing shortwave radiation (Kup), (i–l) incoming longwave radiation (Lup), (m–p) outgoing longwave radiation (Lup), and (q–t) net radiation (QN) from May 2010 to June 2011. The shaded area denotes the standard deviation. Note that while Kdown is input, the rest are model output.


Appendix D: Fitting the maximum photosynthetic rate for the grass and lawn vegetation type

In order to find the maximum photosynthetic rate for the vegetation type of grass to be used in SUEWS simulation, the environmental response functions g(Tair), gq), gθ), and g(Kdown) in Eq. (14) were fitted to observations from an eddy covariance (EC) station (601116.02′′ N, 244956.85′′ E) situated on an urban lawn in Espoo, Finland. A 1.2 m high EC tower was located at the center of the urban lawn covering an area of 0.7 ha. The EC setup consisted of a three-dimensional sonic anemometer (Metek GmbH, Germany) for measuring the three wind components and sonic temperature and a closed-path infrared gas analyzer (LI-7200; LI-COR, Lincoln, NE, USA) for measuring the CO2 and H2O mixing ratios. The gas analyzer inlet was positioned 13 cm below the anemometer, and air was drawn into the gas analyzer using a 60 cm length of steel tube, having an inner diameter of 4.57 mm and a mean flow rate of 12 L min−1. The tube was heated to avoid water vapor condensation on tube walls. The raw EC data were sampled at 10 Hz and stored for post-processing. The steps before 30 min flux calculations consisted of despiking, linear detrending, and planar fitting of the raw data.

The biogenic CO2 flux FC,bio from EC measurements was partitioned into Fres and Fpho using the nighttime temperature dependency of Fres. The nighttime FC,bio was considered nighttime Fres, and it was related to the observed Tair by fitting the exponential model:

(D1) F res = a grass exp ( T air b grass ) ,

where agrass and bgrass are parameters controlling the temperature dependency. The nighttime temperature dependency of Fres was then extrapolated to daytime, and Fpho was then calculated by

(D2) F pho = F C , bio - F res .

Then Fpho was used as a dependent variable, whereas on-site measurements of net radiation (CNR4; Kipp&Zonen, Delft, Netherlands), air temperature and relative humidity (HMP110 A15; Vaisala Oyj, Vantaa, Finland), and soil moisture (ML3 Thetaprobe; Delta-T, Cambridge, UK) were used to estimate the independent variables in g(Tair), gq), gθ), and g(Kdown). Atmospheric pressure from the Finnish Meteorological Institute Kumpula station was used to calculate Δq. Additional reference values of soil properties (field capacity and wilting point), which were estimated to be same as in an urban lawn in Kumpula (Järvi et al.2019), were used to calculate ΔθWP. Parameters Fpho,max, and G2G6 were fitted using a nonlinear least-squares approach.

Data from mid-July to the end of August 2021 were used and in the fitting only data points with Kdown< 10 W m−2 and Δq > 1 g kg−1 were selected (Havu et al.2022a). After applying a bootstrapping method to randomly select seven-eighths of the data 100 times, the final parameters fitted for grass were obtained as medians with the uncertainties as follows: Fpho,max= 5.497 ± 0.110 µmolm-2s-1, G2= 195.019 ± 5.601 W m−2, G3= 0.741 ± 0.008, G4= 0.413 ± 0.015, G5= 30.000 ± 0.000 C, and G6= 0.500 ± 0.000 mm−1. The value of Fpho,max is used for grass and lawn in the Beijing simulation.

Appendix E: Maximum conductance of urban green space in Beijing

The weighted average of gmax is 7.0 mm s−1 for the deciduous tree weighted by the population ratio of each species (Table E1). The average gmax for grass is 3.7 mm s−1. Note that when the species in the modeled area are known, we suggest the gmax is selected accordingly. Here, we seek values that can represent the overall vegetation over the urban area in Beijing, and therefore the average is taken from different species.

Xu et al. (2020)Wang et al. (2018)Xu et al. (2020)Song et al. (2015)Chen et al. (2021)Wang et al. (2006)Wang et al. (2006)Wang et al. (2006)

Table E1Maximum conductance (gmax) for each vegetated surface.

a Obtained from a field survey over Beijing (Ma et al.2019). b obtained by dividing maximum canopy conductance by G1 (= 3.5)

Download Print Version | Download XLSX

Appendix F: Comparison of parameters for modeling domain with different radii

Table F1Land surface fraction for modeling domain with different radii, where “Paved” denotes paved surface, “Bldgs” denotes buildings, “EveTr” denotes evergreen trees, “DecTr” denotes deciduous trees, “Grass” denotes grassland and lawn, “Bsoil” denotes bare soil, and “Water” denotes waterbodies.

Download Print Version | Download XLSX

Figure F1Diurnal cycles of (a) traffic rate and (b) population density for the modeling domain with different radii, where WD denotes weekday and WE denotes weekday.


Code and data availability

The datasets are openly available, including the complete model cases of SUEWS; the meteorological, radiation, and turbulent flux measurements; codes for LAI model optimization using CMA-ES; and codes to reproduce the statistics and figures ( et al.2022).


The supplement related to this article is available online at:

Author contributions

YZ, MH, HL, and LJ conceptualized the study. YZ performed the SUEWS model cases, analyzed the results, and prepared the figures. HL, XC, YW, and JA collected the data. HSL performed the fitting of maximum photosynthetic rate. LJ and HL supervised the study. All authors contributed to writing and preparing the manuscript.

Competing interests

At least one of the (co-)authors is a member of the editorial board of Geoscientific Model Development. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


This study is funded by National Natural Science Foundation of China (grant no. 42161144010), the China Scholarship Council (grant no. 202104910363), the Tiina and Antti Herlin Foundation (grant no. 20200027), the Academy of Finland (grant nos. 321527 and 337549), and the Strategic Research Council established within the Academy of Finland (grant no. 335201). We sincerely appreciate Weiyi Zuo from Institute of Acoustics at the Chinese Academy of Sciences for his significant technical support in realizing CMA-ES method in this study.

Financial support

This research has been supported by the National Natural Science Foundation of China (grant no. 42161144010), the China Scholarship Council (grant no. 202104910363), the Tiina ja Antti Herlinin säätiö (grant no. 20200027), the Academy of Finland (grant nos. 321527 and 337549), and the Academy of Finland, Strategic Research Council (grant no. 335201).

Review statement

This paper was edited by Hisashi Sato and reviewed by two anonymous referees.


Alexander, P. J., Bechtel, B., Chow, W. T., Fealy, R., and Mills, G.: Linking urban climate classification with an urban energy and water budget model: Multi–site and multi–seasonal evaluation, Urban Climate, 17, 196–215,, 2016. a

Alvarez, R. and Weilenmann, M.: Effect of low ambient temperature on fuel consumption and pollutant and CO2 emissions of hybrid electric vehicles in real-world conditions, Fuel, 97, 119–124,, 2012. a

Ao, X., Grimmond, C., Liu, D., Han, Z., Hu, P., Wang, Y., Zhen, X., and Tan, J.: Radiation fluxes in a business district of Shanghai, China, J. Appl. Meteorol. Clim., 55, 2451–2468,, 2016. a, b

Ao, X., Grimmond, C., Ward, H., Gabey, A., Tan, J., Yang, X.-Q., Liu, D., Zhi, X., Liu, H., and Zhang, N.: Evaluation of the Surface Urban Energy and Water balance Scheme (SUEWS) at a dense urban site in Shanghai: Sensitivity to anthropogenic heat and irrigation, J. Hydrometeorol., 19, 1983–2005,, 2018. a, b, c

Awal, M., Ohta, T., Matsumoto, K., Toba, T., Daikoku, K., Hattori, S., Hiyama, T., and Park, H.: Comparing the carbon sequestration capacity of temperate deciduous forests between urban and rural landscapes in central Japan, Urban For. Urban Gree., 9, 261–270,, 2010. a

Best, M. and Grimmond, C.: Key conclusions of the first international urban land surface model comparison project, B. Am. Meteorol. Soc., 96, 805–819,, 2015. a

Best, M. and Grimmond, C.: Modeling the partitioning of turbulent fluxes at urban sites with varying vegetation cover, J. Hydrometeorol., 17, 2537–2553,, 2016. a

Björkegren, A. and Grimmond, C.: Net carbon dioxide emissions from central London, Urban Climate, 23, 131–158,, 2018. a, b

BMBS: Beijing Statistical Yearbook, Beijing Municipal Bureau of Statistics, China Statistics Press, Beijing, ISBN 7-5037-0802-6, 2017. a, b

Boegh, E., Soegaard, H., Broge, N., Hasager, C., Jensen, N., Schelde, K., and Thomsen, A.: Airborne multispectral data for quantifying leaf area index, nitrogen concentration, and photosynthetic efficiency in agriculture, Remote Sens. Environ., 81, 179–193,, 2002. a

Chen, F., Kusaka, H., Bornstein, R., Ching, J., Grimmond, C., Grossman-Clarke, S., Loridan, T., Manning, K. W., Martilli, A., and Miao, S.: The integrated WRF/urban modelling system: development, evaluation, and applications to urban environmental problems, Int. J. Climatol., 31, 273–288,, 2011. a

Chen, S., Chen, Z., and Zhang, Z.: The canopy stomatal conductance characteristics of Pinus tabulaeformis and Acer truncatum and their environmental responses in the mountain area of Beijing, Chinese Journal of Plant Ecology, 45, 1329–1340,, 2021 (in Chinese). a

Cheng, X., Liu, X., Liu, Y., and Hu, F.: Characteristics of CO2 concentration and flux in the Beijing urban area, J. Geophys. Res.-Atmos., 123, 1785–1801,, 2018. a, b, c

Christen, A., Coops, N., Crawford, B., Kellett, R., Liss, K., Olchovski, I., Tooke, T., Van Der Laan, M., and Voogt, J.: Validation of modeled carbon-dioxide emissions from an urban neighborhood with direct eddy-covariance measurements, Atmos. Environ., 45, 6057–6069,, 2011. a, b, c, d

Chu, H., Luo, X., Ouyang, Z., Chan, W. S., Dengel, S., Biraud, S. C., Torn, M. S., Metzger, S., Kumar, J., and Arain, M. A.: Representativeness of Eddy-Covariance flux footprints for areas surrounding AmeriFlux sites, Agr. Forest Meteorol., 301, 108350,, 2021. a

Crawford, B. and Christen, A.: Spatial source attribution of measured urban eddy covariance CO2 fluxes, Theor. Appl. Climatol., 119, 733–755,, 2015. a

Cucchi, M., Weedon, G. P., Amici, A., Bellouin, N., Lange, S., Müller Schmied, H., Hersbach, H., and Buontempo, C.: Near surface meteorological variables from 1979 to 2019 derived from bias-corrected reanalysis, Version 2.0, Copernicus Climate Change Service (C3S) Climate Data Store (CDS) [data set], (last access: 11 May 2022), 2021. a, b

Cui, Y., Zhang, W., Wang, C., Streets, D. G., Xu, Y., Du, M., and Lin, J.: Spatiotemporal dynamics of CO2 emissions from central heating supply in the North China Plain over 2012–2016 due to natural gas usage, Appl. Energ., 241, 245–256,, 2019. a, b, c

Demuzere, M., Harshan, S., Järvi, L., Roth, M., Grimmond, C., Masson, V., Oleson, K., Velasco, E., and Wouters, H.: Impact of urban canopy models and external parameters on the modelled urban energy balance in a tropical city, Q. J. Roy. Meteor. Soc., 143, 1581–1596,, 2017. a

Dou, J., Grimmond, S., Cheng, Z., Miao, S., Feng, D., and Liao, M.: Summertime surface energy balance fluxes at two Beijing sites, Int. J. Climatol., 39, 2793–2810,, 2019. a

Du, M., Wang, X., Peng, C., Shan, Y., Chen, H., Wang, M., and Zhu, Q.: Quantification and scenario analysis of CO2 emissions from the central heating supply system in China from 2006 to 2025, Appl. Energ., 225, 869–875,, 2018. a

Falge, E., Baldocchi, D., Olson, R., Anthoni, P., Aubinet, M., Bernhofer, C., Burba, G., Ceulemans, R., Clement, R., and Dolman, H.: Gap filling strategies for defensible annual sums of net ecosystem exchange, Agr. Forest Meteorol., 107, 43–69,, 2001. a

Fiorella, R. P., Bares, R., Lin, J. C., Ehleringer, J. R., and Bowen, G. J.: Detection and variability of combustion-derived vapor in an urban basin, Atmos. Chem. Phys., 18, 8529–8547,, 2018. a

Foken, T. and Wichura, B.: Tools for quality assessment of surface-based flux measurements, Agr. Forest Meteorol., 78, 83–105,, 1996. a

Fontaras, G., Zacharof, N.-G., and Ciuffo, B.: Fuel consumption and CO2 emissions from passenger cars in Europe–Laboratory versus real-world emissions, Prog. Energ. Combust., 60, 97–131,, 2017. a

Friedl, M. and Sulla-Menashe, D.: MCD12Q1 MODIS/Terra+Aqua Land Cover Type Yearly L3 Global 500 m SIN Grid V006, USGS [data set], (last access: 19 April 2022), 2019. a

Gong, P., Chen, B., Li, X., Liu, H., Wang, J., Bai, Y., Chen, J., Chen, X., Fang, L., and Feng, S.: Mapping essential urban land use categories in China (EULUC-China): Preliminary results for 2018, Sci. Bull., 65, 182–187,, 2020. a

Grimmond, C. and Oke, T.: Heat storage in urban areas: Local-scale observations and evaluation of a simple model, J. Appl. Meteorol., 38, 922–940,<0922:HSIUAL>2.0.CO;2, 1999. a

Grimmond, C. and Oke, T.: Progress in measuring and observing the urban atmosphere, Theor. Appl. Climatol., 84, 3–22,, 2006. a

Grimmond, C. S. B., Blackett, M., Best, M., Barlow, J., Baik, J., Belcher, S., Bohnenstengel, S., Calmet, I., Chen, F., and Dandou, A.: The International Urban Energy Balance Models Comparison Project: First Results from Phase 1, J. Appl. Meteorol. Clim., 49, 1268–1292,, 2010. a

Hansen, N., Auger, A., Ros, R., Finck, S., and Pošík, P.: Comparing results of 31 algorithms from the black-box optimization benchmarking BBOB-2009, in: Proceedings of the 12th annual conference companion on Genetic and evolutionary computation, Association for Computing Machinery, New York, NY, USA, 1689–1696,, 2010. a

Havu, M., Kulmala, L., Kolari, P., Vesala, T., Riikonen, A., and Järvi, L.: Carbon sequestration potential of street tree plantings in Helsinki, Biogeosciences, 19, 2121–2143,, 2022a. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p

Havu, M., Lee, H. S., Soininen, J., and Järvi, L.: Spatial variability of biogenic CO2 flux in Helsinki in 2020 (Version 1), Zenodo [data set],, 2022b. a

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J.-N.: The ERA5 global reanalysis, Q. J. Roy. Meteor. Soc., 146, 1999–2049,, 2020. a

Huete, A.: A comparison of vegetation indices over a global set of TM images for EOS-MODIS, Remote Sens. Environ., 59, 440–451,, 1997. a

Järvi, L., Grimmond, C., and Christen, A.: The Surface Urban Energy and Water Balance Scheme (SUEWS): Evaluation in Los Angeles and Vancouver, J. Hydrol, 411, 219–237,, 2011. a, b, c, d, e, f, g, h, i, j, k

Järvi, L., Grimmond, C. S. B., Taka, M., Nordbo, A., Setälä, H., and Strachan, I. B.: Development of the Surface Urban Energy and Water Balance Scheme (SUEWS) for cold climate cities, Geosci. Model Dev., 7, 1691–1711,, 2014. a, b, c, d

Järvi, L., Grimmond, C., McFadden, J., Christen, A., Strachan, I., Taka, M., Warsta, L., and Heimann, M.: Warming effects on the urban hydrology in cold climate regions, Sci. Rep.-UK, 7, 1–8,, 2017. a

Järvi, L., Havu, M., Ward, H. C., Bellucco, V., McFadden, J. P., Toivonen, T., Heikinheimo, V., Kolari, P., Riikonen, A., and Grimmond, C. S. B.: Spatial modeling of local-scale biogenic and anthropogenic carbon dioxide emissions in Helsinki, J. Geophys. Res.-Atmos., 124, 8363–8384,, 2019. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o

Jiang, X., Zhang, C., Gao, H., and Miao, S.: Impacts of urban albedo change on urban heat island in Beijing-A case study, Acta Meteorol. Sin., 65, 301–307,, 2007 (in Chinese). a

Johansson, E. and Emmanuel, R.: The influence of urban design on outdoor thermal comfort in the hot, humid city of Colombo, Sri Lanka, Int. J. Biometeorol., 51, 119–133,, 2006. a

Karsisto, P., Fortelius, C., Demuzere, M., Grimmond, C. S. B., Oleson, K., Kouznetsov, R., Masson, V., and Järvi, L.: Seasonal surface urban energy balance and wintertime stability simulated using three land-surface models in the high-latitude city Helsinki, Q. J. Roy. Meteor. Soc., 142, 401–417,, 2016. a

Kokkonen, T., Grimmond, C. S. B., Räty, O., Ward, H., Christen, A., Oke, T., Kotthaus, S., and Järvi, L.: Sensitivity of Surface Urban Energy and Water Balance Scheme (SUEWS) to downscaling of reanalysis forcing data, Urban Climate, 23, 36–52,, 2018. a

Kokkonen, T. V., Grimmond, S., Murto, S., Liu, H., Sundström, A.-M., and Järvi, L.: Simulation of the radiative effect of haze on the urban hydrological cycle using reanalysis data in Beijing, Atmos. Chem. Phys., 19, 7001–7017,, 2019. a, b, c, d, e, f, g, h

Konopka, J., Heusinger, J., and Weber, S.: Extensive Urban Green Roof Shows Consistent Annual Net Uptake of Carbon as Documented by 5 Years of Eddy-Covariance Flux Measurements, J. Geophys. Res.-Biogeo., 126, e2020JG005879,, 2021. a

Laine, A., Riutta, T., Juutinen, S., Väliranta, M., and Tuittila, E.-S.: Acknowledging the spatial heterogeneity in modelling/reconstructing carbon dioxide exchange in a northern aapa mire, Ecol. Model., 220, 2646–2655,, 2009. a

Liu, H. Z., Feng, J. W., Järvi, L., and Vesala, T.: Four-year (2006–2009) eddy covariance measurements of CO2 flux over an urban area in Beijing, Atmos. Chem. Phys., 12, 7881–7892,, 2012. a, b, c, d, e, f, g

Liu, S., Pang, H., Zhang, N., Xing, M., Wu, S., and Hou, S.: Temporal variations of the contribution of combustion-derived water vapor to urban humidity during winter in Xi'an, China, Sci. Total Environ., 830, 154711,, 2022. a

Liu, Y., Liu, H., Du, Q., and Xu, L.: Multi-level CO2 fluxes over Beijing megacity with the eddy covariance method, Atmospheric and Oceanic Science Letters, 14, 100079,, 2021. a

Loridan, T., Grimmond, C., Offerle, B. D., Young, D. T., Smith, T. E., Järvi, L., and Lindberg, F.: Local-scale urban meteorological parameterization scheme (LUMPS): longwave radiation parameterization and seasonality-related developments, J. Appl. Meteorol. Clim., 50, 185–202,, 2011. a, b

Lu, P., Yu, Q., Liu, J., and Lee, X.: Advance of tree-flowering dates in response to urban climate change, Agr. Forest Meteorol., 138, 120–131,, 2006. a

Luo, Z., Sun, O. J., Ge, Q., Xu, W., and Zheng, J.: Phenological responses of plants to climate change in an urban environment, Ecol. Res., 22, 507–514,, 2007. a

Ma, J.: The study on urban forest structure and eco-service in the Sixth Ring Road of Beijing, Thesis, Chinese Academy of Forestry, Beijing,, 2019. a

Ma, J., Jia, B., Zhang, W., Liu, X., Li, X., and Liu, J.: The characteristics of urban forest structure within the Sixth Ring Road of Beijing, Chinese Journal of Ecology, 38, 2318–2325,, 2019. a

Marcotullio, P. J., Sarzynski, A., Albrecht, J., Schulz, N., and Garcia, J.: The geography of global urban greenhouse gas emissions: An exploratory analysis, Climatic Change, 121, 621–634,, 2013. a

Masek, J., Vermote, E., Saleous, N., Wolfe, R., Hall, F., Huemmrich, K., Gao, F., Kutler, J., and Lim, T.-K.: A Landsat Surface Reflectance Dataset for North America, 1990–2000, IEEE Geosci. Remote S., 3, 68–72,, 2006. a, b

Masson, V., Le Moigne, P., Martin, E., Faroux, S., Alias, A., Alkama, R., Belamari, S., Barbu, A., Boone, A., Bouyssel, F., Brousseau, P., Brun, E., Calvet, J.-C., Carrer, D., Decharme, B., Delire, C., Donier, S., Essaouini, K., Gibelin, A.-L., Giordani, H., Habets, F., Jidane, M., Kerdraon, G., Kourzeneva, E., Lafaysse, M., Lafont, S., Lebeaupin Brossier, C., Lemonsu, A., Mahfouf, J.-F., Marguinaud, P., Mokhtari, M., Morin, S., Pigeon, G., Salgado, R., Seity, Y., Taillefer, F., Tanguy, G., Tulet, P., Vincendon, B., Vionnet, V., and Voldoire, A.: The SURFEXv7.2 land and ocean surface platform for coupled or offline simulation of earth surface variables and fluxes, Geosci. Model Dev., 6, 929–960,, 2013. a

MHURD: China urban-rural construction statistical yearbook 2016, Ministry of Housing and Urban-Rural Development of the People's Republic of China, China Statistics Press, Beijing, ISBN 9787503779596, 2018. a, b, c

Miao, S., Dou, J., Chen, F., Li, J., and Li, A.: Analysis of observations on the urban surface energy balance in Beijing, Science China Earth Sciences, 55, 1881–1890,, 2012. a

Moriwaki, R. and Kanda, M.: Seasonal and diurnal fluxes of radiation, heat, water vapor, and carbon dioxide over a suburban area, J. Appl. Meteorol., 43, 1700–1710,, 2004. a, b, c

Nordbo, A., Karsisto, P., Matikainen, L., Wood, C. R., and Järvi, L.: Urban surface cover determined with airborne lidar at 2 m resolution–implications for surface energy balance modelling, Urban Climate, 13, 52–72,, 2015. a

Offerle, B., Grimmond, C., and Oke, T. R.: Parameterization of net all-wave radiation for urban areas, J. Appl. Meteorol. Clim., 42, 1157–1173,<1157:ponarf>;2, 2003. a, b

Oke, T.: The Heat Island of the Urban Boundary Layer: Characteristics, Causes and Effects, Springer, Dordrecht,, pp. 81–107, 1995. a

Omidvar, H., Sun, T., Grimmond, S., Bilesbach, D., Black, A., Chen, J., Duan, Z., Gao, Z., Iwata, H., and McFadden, J. P.: Surface Urban Energy and Water Balance Scheme (v2020a) in vegetated areas: parameter derivation and performance evaluation using FLUXNET2015 dataset, Geosci. Model Dev., 15, 3041–3078,, 2022. a, b, c, d

Pataki, D. E., Carreiro, M. M., Cherrier, J., Grulke, N. E., Jennings, V., Pincetl, S., Pouyat, R. V., Whitlow, T. H., and Zipperer, W. C.: Coupling biogeochemical cycles in urban environments: ecosystem services, green solutions, and misconceptions, Front. Ecol. Environ., 9, 27–36,, 2011. a

Peters, E. B. and McFadden, J. P.: Continuous measurements of net CO2 exchange by vegetation and soils in a suburban landscape, J. Geophys. Res.-Biogeo., 117, G03005,, 2012. a

Prairie, Y. T. and Duarte, C. M.: Direct and indirect metabolic CO2 release by humanity, Biogeosciences, 4, 215–217,, 2007. a

Qian, Y., Chakraborty, T., Li, J., Li, D., He, C., Sarangi, C., Chen, F., Yang, X., and Leung, L. R.: Urbanization Impact on Regional Climate and Extreme Weather: Current Understanding, Uncertainties, and Future Research Directions, Adv. Atmos. Sci., 39, 819–860,, 2022. a

Rebmann, C., Göckede, M., Foken, T., Aubinet, M., Aurela, M., Berbigier, P., Bernhofer, C., Buchmann, N., Carrara, A., and Cescatti, A.: Quality analysis applied on eddy covariance measurements at complex forest sites using footprint modelling, Theor. Appl. Climatol., 80, 121–141,, 2005. a

Salmon, O. E., Shepson, P. B., Ren, X., Marquardt Collow, A. B., Miller, M. A., Carlton, A. G., Cambaliza, M. O., Heimburger, A., Morgan, K. L., and Fuentes, J. D.: Urban emissions of water vapor in winter, J. Geophys. Res.-Atmos., 122, 9467–9484,, 2017. a

Sarangi, C., Tripathi, S., Qian, Y., Kumar, S., and Ruby Leung, L.: Aerosol and urban land use effect on rainfall around cities in Indo–Gangetic Basin from observations and cloud resolving model simulations, J. Geophys. Res.-Atmos., 123, 3645–3667,, 2018. a

Song, Y., Li, F., Wang, X., Xu, C., Zhang, J., Liu, X., and Zhang, H.: The effects of urban impervious surfaces on eco-physiological characteristics of Ginkgo biloba: A case study from Beijing, China, Urban For. Urban Gree., 14, 1102–1109,, 2015. a

Stagakis, S., Feigenwinter, C., Vogt, R., and Kalberer, M.: A high-resolution monitoring approach of urban CO2 fluxes. Part 1 – bottom-up model development, Sci. Total Environ., 160216,, 2022. a, b

Sun, T. and Grimmond, S.: A Python-enhanced urban land surface model SuPy (SUEWS in Python, v2019.2): development, deployment and demonstration, Geosci. Model Dev., 12, 2781–2795,, 2019. a

Tan, J., Zheng, Y., Tang, X., Guo, C., Li, L., Song, G., Zhen, X., Yuan, D., Kalkstein, A. J., and Li, F.: The urban heat island and its impact on heat waves and human health in Shanghai, Int. J. Biometeorol., 54, 75–84,, 2010. a

United Nations Department of Economic and Social Affairs: World Urbanization Prospects: The 2018 Revision, United Nations,, 2019. a

Velasco, E. and Roth, M.: Cities as net sources of CO2: Review of atmospheric CO2 exchange in urban environments measured by eddy covariance technique, Geography Compass, 4, 1238–1259,, 2010. a

Wang, L., Fan, S., Hu, F., Miao, S., Yang, A., Li, Y., Liu, J., Liu, C., Chen, S., Ho, H. C., Duan, Z., Gao, Z., and Yang, Y.: Vertical gradient variations in radiation budget and heat fluxes in the urban boundary layer: A comparison study between polluted and clean air episodes in Beijing during winter, J. Geophys. Res.-Atmos., 125, e2020JD032478,, 2020. a

Wang, Y., Yang, J., and Han, L.: Effects of different irrigation volume on photosynthesis of turfgrass (in Chinese), Journal of Beijing Forestry University, 28, 26–31, 2006. a, b, c

Wang, Y., Li, G., Di, N., Clothier, B., Duan, J., Li, D., Jia, L., Xi, B., and Ma, F.: Leaf phenology variation within the canopy and its relationship with the transpiration of Populus tomentosa under plantation conditions, Forests, 9, 603,, 2018. a

Wang, Y., Chang, Q., and Li, X.: Promoting sustainable carbon sequestration of plants in urban greenspace by planting design: A case study in parks of Beijing, Urban For. Urban Gree., 64, 127291,, 2021. a

Ward, H. C., Evans, J. G., and Grimmond, C. S. B.: Multi-season eddy covariance observations of energy, water and carbon fluxes over a suburban area in Swindon, UK, Atmos. Chem. Phys., 13, 4645–4666,, 2013. a

Ward, H. C., Kotthaus, S., Järvi, L., and Grimmond, C. S. B.: Surface Urban Energy and Water Balance Scheme (SUEWS): development and evaluation at two UK sites, Urban Climate, 18, 1–32,, 2016. a, b, c, d, e, f, g, h, i, j, k, l, m

Watts, N., Adger, W. N., Agnolucci, P., Blackstock, J., Byass, P., Cai, W., Chaytor, S., Colbourn, T., Collins, M., and Cooper, A.: Health and climate change: policy responses to protect public health, Lancet, 386, 1861–1914,, 2015. a

Wen, Y., Zhang, S., Zhang, J., Bao, S., Wu, X., Yang, D., and Wu, Y.: Mapping dynamic road emissions for a megacity by using open-access traffic congestion index data, Appl. Energ., 260, 114357,, 2020. a, b

Wen, Y., Wu, R., Zhou, Z., Zhang, S., Yang, S., Wallington, T. J., Shen, W., Tan, Q., Deng, Y., and Wu, Y.: A data-driven method of traffic emissions mapping with land use random forest models, Appl. Energ., 305, 117916,, 2022. a, b

Xu, Y., Li, S., Yuan, X., and Feng, Z.: Emission characteristics of biogenic volatile compounds (BVOCs) from common greening tree species in northern China and their correlations with photosynthetic parameters, Environmental Science, 41, 3518–3526,, 2020 (in Chinese). a, b

Yang, D., Zhang, S., Niu, T., Wang, Y., Xu, H., Zhang, K. M., and Wu, Y.: High-resolution mapping of vehicle emissions of atmospheric pollutants based on large-scale, real-world traffic datasets, Atmos. Chem. Phys., 19, 8831–8843,, 2019. a, b, c

Zhang, D., Gao, J., Tang, D., Wu, X., Shi, J., Chen, J., Peng, Y., Zhang, S., and Wu, Y.: Switching on auxiliary devices in vehicular fuel efficiency tests can help cut CO2 emissions by millions of tons, One Earth, 4, 135–145,, 2021. a

Zhang, H., Zhou, L., Huang, X., and Zhang, X.: Decarbonizing a large City's heating system using heat pumps: A case study of Beijing, Energy, 186, 115820,, 2019. a, b

Zhang, S., Wu, Y., Liu, H., Huang, R., Yang, L., Li, Z., Fu, L., and Hao, J.: Real-world fuel consumption and CO2 emissions of urban public buses in Beijing, Appl. Energ., 113, 1645–1655,, 2014. a, b

Zhang, X., Mi, F., Lu, N., Yan, N., Kuglerova, L., Yuan, S., Peng, Q., and Ma, O. Z.: Green space water use and its impact on water resources in the capital region of China, Phys. Chem. Earth Pts. A/B/C, 101, 185–194,, 2017.  a

Zhang, Y., Yin, P., Li, X., Niu, Q., Wang, Y., Cao, W., Huang, J., Chen, H., Yao, X., and Yu, L.: The divergent response of vegetation phenology to urbanization: A case study of Beijing city, China, Sci. Total Environ., 803, 150079,, 2022. a

Zhao, X., Zhou, Y., Chen, W., Li, X., Li, X., and Li, D.: Mapping hourly population dynamics using remotely sensed and geospatial data: a case study in Beijing, China, GISci. Remote Sens., 58, 717–732,, 2021. a

Zheng, Y., Liu, H., and Cheng, X.: Datasets for simulating heat and CO2 fluxes in Beijing using SUEWS V2020b, Zenodo [data set],, 2022. a, b, c

Zhu, X., Ni, G., Cong, Z., Sun, T., and Li, D.: Impacts of surface heterogeneity on dry planetary boundary layers in an urban-rural setting, J. Geophys. Res.-Atmos., 121, 12164–12179,, 2016. a

Short summary
The performance of the Surface Urban Energy and Water Balance Scheme (SUEWS) is evaluated against the observed surface exchanges (fluxes) of heat and carbon dioxide in a densely built neighborhood in Beijing. The heat flux modeling is noticeably improved by using the observed maximum conductance and by optimizing the vegetation phenology modeling. SUEWS also performs well in simulating carbon dioxide flux.