Simulating heat and CO 2 fluxes in Beijing using SUEWS V2020b: Sensitivity to vegetation phenology and maximum conductance

. The Surface Urban Energy and Water Balance Scheme (SUEWS) has recently been introduced to include a bottom-up approach to modelling carbon dioxide (CO 2 ) emissions and sink in urban areas. In this study, SUEWS is evaluated against radiation flux observations and eddy covariance (EC) measured turbulent fluxes of sensible heat ( Q H ), latent heat ( Q E ), and CO 2 ( F C ) at a densely built neighborhood in Beijing. The model sensitivity to maximum conductance ( g max ) and leaf area index (LAI) is examined. Site-specific g max is obtained from observations over local vegetation species, and LAI parameters 5 by optimization with remotely sensed LAI obtained from a MODIS/Terra data product. For simulation of anthropogenic CO 2 components, local traffic and population data are collected. In model evaluation, the mismatch between the measurement source area and simulation domain is also considered. Using the optimized g max and LAI, the modelling of heat fluxes is noticeably improved, showing higher correlation with observations, lower bias, and more realistic seasonal dynamics of Q E and Q H . In comparison to heat fluxes, the F C module 10 shows lower sensitivity to the choice of g max and LAI. This can be explained by the low relative contribution of vegetation to net F C in the modelled area. SUEWS successfully reproduces the average diurnal cycle of F C and annual cumulative sums. Depending on the size of the simulation domain, the modelled annual accumulated F C ranges from 7.2 to 8.5 kg C m − 2 yr − 1 , when compared to 7.5 kg C m − 2 yr − 1 observed by EC. Traffic is the dominant CO 2 source, contributing 63–73 % to the annual total CO 2 emissions, followed by human metabolism (14–18%), respiration released by vegetation and soil (6– 15 11%) and building heating (6–9%). Vegetation photosynthesis offsets only 4–8 % of the total CO 2 emissions. We highlight the importance of choosing optimal LAI parameters and g max when SUEWS is used to model surface fluxes. The F C module of SUEWS is a promising tool in quantifying urban CO 2 emissions at the local scale, and therefore assisting to mitigate urban CO 2 emissions.


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 21 st century (United Nations Department of Economic and Social Affairs, 2019).The urban expansion reshapes the morphological, thermal, and dynamical properties of the land surface (Grimmond and Oke, 2006;Oke, 1995;Zhu et al., 2016), and intensive human activities contribute to a major source of greenhouse gas emissions (Marcotullio et al., 2013;Velasco and Roth, 2010).This has influenced urban climate from micro to regional scales (Johansson and Emmanuel, 2006;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-atmosphere interactions, including the exchanges of energy, water and CO 2 , and hydrological processes (Chen et al., 2011;Masson et al., 2013).The results from the First International Urban Land Surface Model Comparison Project suggest that the most important processes for surface energy balance in the urban environment were radiative and vegetation processes (e.g.surface fraction of vegetation, seasonal cycle of vegetation phenology) (Grimmond et al., 2010;Best and Grimmond, 2015;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 sensitive in such environments (Best and Grimmond, 2016).
The Surface Urban Energy and Water balance Scheme (SUEWS) is one of the widely tested ULSMs (Järvi et al., 2011(Järvi et al., , 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 Grimmond, 2019).SUEWS has demonstrated good performance against hydrological observations and surface flux observations in several cities in Europe, North America, and Asia (Järvi et al., 2011;Alexander et al., 2016;Ward et al., 2016;Ao et al., 2018;Havu et al., 2022a).In SUEWS, 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 modelled LAI leads to over-or underestimation in Q E or Q H (Ao et al., 2018;Ward et al., 2016).They highlight 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 outskirt of cities. Apart from LAI, the maximum conductance (g max ) is also critical in scaling the surface conductance (g s ), and therefore the available energy distribution (Ward et al., 2016).However, the impact of LAI-related parameters and g max on modelled turbulent fluxes has received insufficient attention in urban areas.
Recently, the module of local-scale CO 2 flux (F C ) was incorporated into SUEWS (Järvi et al., 2019).It was found to give reasonable annual sum, seasonal and diurnal cycles against observed F C in Helsinki, Finland, suggesting that the bottomup CO 2 emission or sink 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 CO 2 emission, quanti-https://doi.org/10.5194/gmd-2022-305Preprint.Discussion started: 9 February 2023 c Author(s) 2023.CC BY 4.0 License.
where LAI max,i and LAI min,i for each vegetation type can be obtained from literature or determined from observations, ω 1/2,GDD/SDD,i represent the growing or senescence rates derived for each study site or use their default values (Järvi et al., 2011;Omidvar et al., 2022), and T d−1 is the previous day air temperature mean.

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

Turbulent heat fluxes
Latent heat flux (Q E , W m −2 ) is calculated using the modified Penman-Monteith equation for each surface type: where Q N (W m −2 ) is the net all-wave radiation, Q F (W m −2 ) the anthropogenic heat flux, ∆Q S (W m −2 ) the net storage heat flux, ρ (kg m −3 ) the air density, c p (J kg −1 K −1 ) the specific heat capacity of air at constant pressure, V P D (Pa) the vapour pressure deficit, s (Pa • C −1 ) the slope of the saturation vapour pressure curve, γ (Pa • C −1 ) the psychrometric constant, r av (s mm −1 ) the aerodynamic resistance for water vapour, and r s (s mm −1 ) the surface resistance.r s , or its inverse surface conductance g s (mm s −1 ), has the form: where g max,i is the maximum conductance of vegetation type i, f r i is the surface fraction of i, G 1 is a constant connecting stomatal conductance to canopy conductance, g(K down ), g(∆q), g(T air ), and g(∆θ) are environmental response functions on K down , specific humidity deficit (∆q), air temperature (T air ), and soil moisture deficit (∆θ), respectively.The functions have forms (Ward et al., 2016): where and Coefficients G 2 − G 6 determine the shape of the curves describing responses of stomatal conductance to each environmental variable.K down,max (W m −2 ) is the maximum incoming solar radiation, T L and T H ( • C) are the lower and upper limits for temperature at which photosynthesis and transpiration are off, and ) is calculated from model input RH, T air ( • C) is either model input or simulated at 2 m height, and ∆θ (mm) is simulated within SUEWS (Järvi et al., 2017).Q H is determined as the residual from the surface energy balance equation:

CO 2 flux
The F C module adopts a bottom-up approach to determine the local-scale F C (µmol m −2 s −1 ), accounting for both anthropogenic (F C,ant ) and biogenic (F C,bio ) components (Järvi et al., 2019): In the Eq. 10, F M are CO 2 emissions from human metabolism, F V emissions from traffic, F B emissions from buildings (e.g., combustion of natural gas, coal and wood), F P emissions from local-scale point sources, F pho photosynthesis, and F C,ant relates to energy balance through Q F .F M and F V are estimated based on an inventory approach, i.e., based on population density or traffic rate, and their emission factors (EFs).Hourly CO 2 emissions from human metabolism on workdays/weekends (F M,h,d , µmol m −2 s −1 ) are calculated using: where p h,d is the population density (cap ha −1 ), H a,h,d activity by hour calculated from the diurnal profiles for population and activity, and C M CO 2 released per person (µmol CO 2 s −1 cap −1 ).C M varies between nighttime minimum and daytime maximum values (C M (min,max) ).
Hourly traffic CO 2 emissions (F V,h,d ) on weekdays or weekends are calculated from Hourly building CO 2 emissions (F B,h,d ) on weekdays or weekends are calculated from where f r heat is the fraction of fossil fuels used for heating, Q F,heat building heat emission at local scale estimated from the heating-degree-day model (Järvi et al., 2011), f r nonheat fraction of fossil fuels used for energy other than heating (e.g.combustion from domestic cooking), Q F,base non-temperature related anthropogenic heat flux (W m −2 ), f r QF,base,BEU,d the fraction of the Q F,base coming from building energy use on weekdays or weekends, and E CO2perJ the EF for fuels in building energy use.Emissions from single-point sources such as power plants and industrial activities can be included as model inputs.
F pho relates to energy balance through LAI and the environmental responses of surface conductance (Eq.3).F pho is calculated using where F pho,max,i is the maximum photosynthetic rate for vegetation type i.
Soil and vegetation respiration F res (µmol m −2 s −1 ) follows an exponential dependency on T air : where a i and b i are parameters controlling the temperature dependency, and 0.6 µmol m −2 s −1 is the minimum respiration in winter time.

Study site and measurements
The model domain is located around a 325-m meteorological tower constructed by Institute of Atmospheric Physics, Chinese Academy of Sciences (IAP tower, 39 • 58' N, 116 • 22' E, 60 m above sea level) within the 6 th Ring area of Beijing, China (Fig. 1 a).An EC setup at the height of 47 m on the IAP tower continuously measures the surface fluxes of Q E , Q H and F C using a 3-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 the IAP tower roughly covers 80% accumulated flux footprint area (Liu et al., 2012).This area is mainly covered by impervious surfaces (Fig. 1 b).Three patches of urban green spaces are situated to the east, south and west to the IAP tower, while the other vegetation is scarcely located along the roads and near the buildings.Most of the impervious surfaces in the source area are residential buildings (Fig. 1 c).A more detailed description of the surroundings and the used instruments can be found in Liu et al. (2012), Liu et al. (2021) and Cheng et al. (2018).with building heights over 50 m (112-128, 160-243, 314-3 • ) are removed (Kokkonen et al., 2019).( 4) Stationarity test: data points with stationary indicator > 30% are filtered out (Foken and Wichura, 1996).The percentages of data removed through these 4 steps are 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 Q E , 7338 for Q H and 7797 for F C .These 30-minute flux observations are resampled to one-hour resolution.In addition, F C is gap-filled using seasonal mean diurnal cycle in order to calculate the seasonal and annual sums (Falge et al., 2001).
4 Model run

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 modelling.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 and at hourly temporal resolution.WFDE5 is evaluated against observed meteorological observations before it is used as the forcing data for SUEWS (Appendix B).

Land cover
Land cover types and their fractions needed in the model runs are estimated based on aerial images.Paved surface accounts 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 6 th Ring area, the population of deciduous species accounts for 82% of the total number of woody plants investigated (Ma, 2019).Therefore, the surface fraction of deciduous trees is set to 11% and evergreen trees 2%.
To calculate the storage heat flux, Objective Hysteresis Model (OHM) is used (Grimmond and Oke, 1999), and the coefficients for paved surface is set to the weighted average of case AN99 to represent the asphalt surface following Ward et al. (2016).

Human activity
In this study, local parameters of traffic, population dynamics, and building energy use are incorporated when compared to Kokkonen et al. (2019) as CO 2 emissions are more sensitive to them when compared to anthropogenic heat flux.Hourly traffic rate on weekdays or weekends is calculated as the sum of traffic volume weighted by road length within the study domain, and mean daily traffic rates (T r wd/we ) and their diurnal profiles are then obtained (Table 1, Fig. 2).The average traffic CO 2 EFs are calculated using the method following Zhang et al. (2018) based on the self-reported fuel records, on-road measurements, vehicle category and speed obtained in Beijing (Zhang et al., 2014;Wen et al., 2020Wen et al., , 2022) ) (Table 1).Traffic rate during daytime is higher than in the nighttime, and it is higher on weekdays than on weekends throughout the day (Fig. 2).
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 years 2015 and 2016 (Zhao et al., 2021) (Fig. 2).As there are several residential building areas located within the model domain (Fig. 1 c), population density increases in the evening when residents get home from work, remains high throughout the night, and decreases in the morning.
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.In 2015, the ratio of district boilers heating capacity to cogeneration plants was 4.2:1 (Zhang et al., 2019).Cogeneration plants are usually located at sub-urban or rural areas, and there are no cogeneration plants within the model domain, so their contribution to CO 2 emission is neglected in this study.In comparison, over 5000 coal-fired and 1000 gas-fired heating boilers are located surrounding the populated areas in 2014 (Cui et al., 2019).The exact number or heating capacity is unknown in the study area, but three boiler plants are seen within 2 km distance from the IAP tower through Baidu Map.Their CO 2 emissions are very likely to be observed by EC during heating season.Therefore, fraction of fossil fuels used for heating (F r heat ) is set to 0.81 based on the ratio of heating capacity.SUEWS first estimates the anthropogenic heat, then converts the heat into CO 2 combustion using the EF and emission fractions.The value of natural gas consumption over the annual total heating supply is 3.2×10 6 ton coal equivalent (tce) and the value for coal consumption 2.6×10 6 tce in 2015.The consumption of rest of the fuel types is only 3.8×10 5 tce (MHURD, 2018).The EFs of heating supply are 96.51Tg CO 2 /10 19 J for coal-fired boiler, and 56.17 Tg CO 2 /10 19 J for gas-fired boiler, respectively (Du et al., 2018).Therefore, the EF for fuels used in building energy use (E CO2perJ ) in SUEWS takes the average of the EFs of natural gas and coal weighted by their annual consumption, i.e., 0.1688 µmol CO 2 J −1 (Table 1).In addition, SUEWS needs a temperature limit (base temperature, T base_HC ) to indicate when heating takes place in the heating-degree-day model.Central heating usually starts around 15 th of November and lasts until 15 th of March, when the outdoor air temperature is around 12 • C. Therefore, this value is given to SUEWS T base_HC to replace the default value (18.2 • C) (Järvi et al., 2011).
Household fuel combustion, mainly from cooking, takes place throughout the year, but its CO 2 emission are relatively small.
An observational study shows that CO 2 emitted from fuel combustion in household cooking contributes only 6% of the indoor CO 2 , and this percentage is low, compared with contribution by human metabolism (30%) (Shen et al., 2020).Therefore, we set 20% to non-heating energy fraction (F r nonheat ).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 modelled area between 2010 and 2016 (figure not shown).Therefore, we assume that the evaluation of radiation fluxes using observations in years 2011-2012 holds true in 2016.

Sensitivity to vegetation-related parameters
In order to test the model's sensitivity of radiation and turbulent fluxes to vegetation-related parameters, four model runs are designed as follows: 1. Case base: control run where model parameters are considered "default" following Kokkonen et al. (2019)     the flux footprint of roughly 60%, 70%, and 80-90%, respectively (Liu et al., 2012).Used vegetation parameters are as in case gs_LAI, but land surface fractions, population, and traffic parameters are modified accordingly (Appendix E).

Statistical metrics for model evaluation
Common statistical metrics are used to quantify the model performance, including, coefficient of determination (R 2 ), rootmean-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 R 2 .The other statistical metrics are defined as follows: where ŷi is the modelled and y i the measured value.Statistical metrics are calculated at annual and seasonal scale, i.e., DJF (winter), MAM (spring), JJA (summer), SON (autumn).

Seasonal dynamics of optimized LAI
The case base has the ability to simulate the onset of leaf growth and the ending of senescence well, but the general vegetation phenology modelling performs poorly (Fig. 3).The performance of LAI modelling is remarkably good after the optimization, Although previous studies suggest that LAI is generally modelled well using default parameters following Järvi et al. (2011), Ward et al. (2016) reported that leaf-on is reached too early and suddenly in spring in two UK cities.In the contrary, Ao et al.
(2018) showed that the simulated LAI might be lower than reality when vegetation remains green in winter and spring in Shanghai, China.These lead to the bias in g s and therefore Q E estimates (Eq.2-3).In Beijing, the rainy season lasts from May to October, while the other time of the year is the dry season (Liu et al., 2012).This influences the LAI behaviour such as the growth and senescence rates due to the lack of soil moisture in spring and autumn (Omidvar et al., 2022).Therefore, it is necessary to evaluate the LAI model when applied to a different city, and optimization of the LAI model is recommended when the LAI model performs poorly.

Evaluation of radiation fluxes
Four model experiments are conducted to examine the impact of g max and LAI seasonal dynamics on the modelled radiation fluxes (See Sect.4.4.1).However, the choice of g max and LAI parameters has only a minor impact on the modelled radiation fluxes (K down , K up , L down , L up ).For all of them, the R 2 between the case gs_LAI and the other cases is 1, and the absolute The model performance of SUEWS in simulating radiation fluxes is good (Fig. 4, Table 5) and it can reproduce the diurnal 300 cycle of each radiation flux well (Fig. 5).K up is overestimated in all seasons with MBE ranging from 0 to 10 W m −2 .At least partly, this overestimation can be explained by the positive bias (MBE > 15 W m −2 in all seasons) of WFDE5 K down when compared to the observed K down (Appendix B).The overestimation in K up 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 average albedo for modelling domain given to SUEWS is 0.14, which is relatively high but still consistent with the observations.Larger positive bias is observed in summer than in winter.In SUEWS, albedo is allowed to vary between summer and winter if specified, and this can help to lower the bias in summer when the albedo for vegetation is low, but this is likely to have only a minor impact because the vegetated fraction is low.We conclude that SUEWS is applicable to provide realistic estimate of radiation fluxes in Beijing, in general accordance 315 with previous studies (Järvi et al., 2014;Karsisto et al., 2016;Ward et al., 2016), despite the absence of site-specific parameters.

Evaluation of turbulent heat fluxes modelling
Both observed Q E and Q H reach their maxima around noon (Fig. 6).The observed Q E has the largest amplitude during the summer months (JJA), while Q H during the spring months (MAM).All four model runs capture their diurnal cycles, but large differences in the amplitude and model performance are observed among the model runs (Fig. 7).
In the case base, the model overestimates Q E (with MBE from -7.4 to 48.6 W m −2 ) except in winter months (DJF) (Fig. 7).
With the optimized LAI (case LAI), slightly better model performance in Q E is seen with lower RMSE (11.5-91.2W m −2 ) and larger R 2 (0.20-0.58) when compared to case base (with RMSE 11.7-96.1 W m −2 and R 2 0.20-0.51).(Fig. 7 a-c).Q E is to a large extent determined by surface conductance which for the modelled area is scaled by g max of each vegetated surface (Eq.3).Clear improvement is observed when local g max is used (case gs), especially during summer, when the overestimation of Q E is largely reduced, RMSE drops to 11.3-54.7 W m −2 , and R 2 increases to 0.25-0.61.The introduction of local g max improves the model performance more than the optimized LAI.With both optimized LAI and local g max introduced (case gs_LAI), the R 2 and RMSE are similar to those at case gs, while MBE drops below 0 in each season, indicating the Q E is slightly underestimated throughout the year.The underestimation can be explained by the missing of combustion-derived water vapor in SUEWS as observational studies have shown that combustion-derived water vapor often contributes 5-10% of total urban humidity during heating season (Fiorella et al., 2018;Liu et al., 2022;Salmon et al., 2017).The observed Q E 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 Q E , might lead to Q E at this magnitude.In SUEWS, Q H is estimated as the residual of energy balance, and is therefore directly affected by the modelled Q E .As a result of overestimating Q E in case base, Q H is greatly underestimated with MBE -51.1-14.4W m −2 and RMSE 60.9-116.0W m −2 .R 2 in summer months and autumn months (SON) are lower than 0.1.The model performance is barely improved by optimized LAI (case LAI), but is noticeably improved after local g max is introduced (case gs) (Fig. 7 d-f).The best model performance for Q H is obtained by case gs_LAI, decreasing the magnitudes of MBE to -7.6-8.9W m −2 and RMSE to 58.0-69.1 W m −2 , and increasing R 2 to 0.42-0.72 (Fig. 7 d-e).
Our results suggest that both altering the values of g max based on observations and LAI model optimization can serve as effective approaches to improving heat fluxes modelling.Expectedly, SUEWS has difficulty in capturing the hourly variability of F C , resulting in the overall low R 2 (0.16-0.21) and high RMSE (12.6-16.4µmol m −2 s −1 ) (Fig. 10 a).On one hand, observed F C has great variability at hourly scale throughout the year indicated by the large deviation (Fig. 8).On the other, under the current parameterization, two of the anthropogenic F C components are static: modelled traffic emission diurnal cycle is only dependent on whether it is weekend or weekday, and modelled human metabolism diurnal cycle is invariable throughout the year (Fig. 9), making it difficult to capture the extreme values.Other urban F C bottom-up modelling studies also reported similar challenges in modelling F C hourly variability in Basel and in Helsinki (Stagakis et al., 2022;Järvi et al., 2019).

Evaluation of CO
SUEWS gives a reasonable estimate of annual accumulated F C (8.4 kg C m −2 yr −1 ), which is 12% higher than the observed gap-filled value (7.5 kg C m −2 yr −1 ).F C is overall slightly underestimated in winter months (with MBE around -0.7 µmol m −2 s −1 ) while overestimated in the other seasons, with the largest MBE (4.6-5.5 µmol m −2 s −1 ) in summer (Fig. 10 c).
SUEWS overestimates summer F C , which might be explained by the underestimation of vegetation photosynthetic rate due to the lack of site-specific parameters.

Model uncertainties
Uncertainties in traffic emission originate from the traffic rates and EFs.SUEWS adopts static traffic EFs, and neglects the relationship between traffic emission and T air as reported by Alvarez and Weilenmann (2012) and Fontaras et al. (2017).In order to examine the impact of seasonal variation of T air 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, https://doi.org/10.5194/gmd-2022-305Preprint.Discussion started: 9 February 2023 c Author(s) 2023.CC BY 4.0 License.
of the current values (Table 1), the human metabolic emission will increase and the annual F C will be 4% higher than the original estimate.
Building emission are calculated based on the Q F estimates and heating fraction.Modelled average Q F in December is 52.7 W m −2 , which is higher than another model estimate (21.6 W m −2 ) in the modelled area (Wang et al., 2020).Observations of Q F are rarely available, and thus these Q F 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 modelled area is unknown.However, the building emission estimate (0.6 kg C m −2 yr −1 ) falls in the range of estimates (~0-3.0kg C m −2 yr −1 ) by other cities (Björkegren and Grimmond, 2018;Järvi et al., 2019;Moriwaki and Kanda, 2004;Christen et al., 2011).
The modelled respiration is larger than CO 2 assimilated through photosynthesis in our study.At the annual scale, urban vegetative surfaces can serve as a net CO 2 sink (Awal et al., 2010;Konopka et al., 2021), but may also vary from a net CO 2 sink to a source (Peters and McFadden, 2012).Admittedly, bias might exist in biogenic CO 2 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 CO 2 sequestrated by the local vegetation, and overestimate net F C in the growing season as shown in Sect.5.4.1.

The impact of the modelling domain size
The surroundings of the IAP tower are heterogeneous in terms of land surface fraction and mean daily traffic rate (Fig. 11 a).
The fraction of vegetated surfaces is higher closer to the tower than further away due to the green spaces adjoining the IAP tower (Fig. 1 b).Additionally, there is a traffic hot spot on the North 3 rd Ring Road located 850 m to the south of IAP tower, 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 modelling domain is 1000 m or 1500 m when compared to domains with lower radii (Fig. 11 a).Thus, the modelled annual accumulated F C largely depends on the modelling domain size chosen, giving estimates of 7.2, 7.4, 8.4 and 8.5 kg C m −2 yr −1 for the radii of 500 m, 750 m, 1000 m, and 1500 m, respectively.
Observational annual F C (7.5 kg C m −2 yr −1 ) falls within this range, which indicates the good model performance of SUEWS (Fig. 11 b).
Land surface fractions are critical parameters in turbulent flux modelling, but the use of site-specific fractions covering a certain area might not yield the best model performance (Demuzere et al., 2017;Loridan and Grimmond, 2012).Approximating the source area by choosing the radius of ≥80% footprint fetch (i.e.≥1000 m in this study) does not give the closest estimate of annual F C , showing agreement with previous studies.The model performance on heat and F C fluxes is reasonably good for each radius (figure not shown), but the selection of radius causes a difference up to 17.5% on annual F C estimate (Fig. 11b).This can be explained by the mismatch between the modelling domain and the real flux source area-the single fixedextent modelled area cannot perfectly represent the land surface characteristics (e.g., the nonuniform land cover and human activities), biasing turbulent flux modelling (Chu et al., 2021;Laine et al., 2009).First, the accumulated footprint area of the observed fluxes is irregular in shape and vary with time (Liu et al., 2012), while many studies (Demuzere et al., 2017 et al., 2019) including this study select a circular area to be modelled and evaluated.Changing the size of modelling domain is challenging when soil processes are taken into account in SUEWS.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).When the modelling 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.We conclude that the land surface model single-point evaluation needs to be performed with the awareness of the mismatch between flux source area and modelling domain, especially over heterogeneous surfaces such as urban surfaces.Regardless of the modelling domain size, traffic is the dominant CO 2 source contributing 63-73% to the total CO 2 emissions, followed by human metabolism (14-18%), respiration released by vegetation and soil (6-11%), and heating (6-9%).Vegetation photosynthesis offsets only 4-8% of the total annual CO 2 emissions (Fig. 12).Several bottom-up modelling studies show that on-road traffic is the greatest source in a densely built neighbourhood, contributing to 70% in central London, 61% in Helsinki, 53%-78% in Tokyo, 70% in Vancouver, while human metabolism also plays an important role, contributing 5-39% the annual total F C (Björkegren and Grimmond, 2018;Järvi et al., 2019;Moriwaki and Kanda, 2004;Christen et al., 2011).Our results are in general agreement with these studies.The contribution of building emission is more variable among cities: a contribution of 70% is reported in Basel (Stagakis et al., 2022), while ~0% in Helsinki (Järvi et al., 2019), and our estimate falls in this range.The direct CO 2 sequestration by urban plants is minor compared with the total CO 2 emissions in this densely built  radii of circular area ranging from 500 m to 1500 m (i.e., accumulated footprint area from 60% to 80-90%), the modelled annual F C ranges from 7.2 to 8.5 kg C m −2 yr −1 , which is comparable with the EC observations (7.5 kg C m −2 yr −1 ).This shows the model performs well also on the annual scale.The variation in annual cumulative F C with radii can be explained by spatial heterogeneity over the urban landscape, and the mismatch between the modelling domain and flux source area.Regardless of the modelling domain size, traffic is the dominant CO 2 source, contributing 63-73% to the total CO 2 emissions, followed by metabolism (14-18%), respiration released by vegetation and soil (6-11%), and heating (6-9%).Vegetation photosynthesis offsets only 4-8% of the CO 2 emissions.
We highlight the importance in choosing more site-specific LAI parameters or g max when using SUEWS in modelling heat fluxes, before more advanced application regarding urban climate and hydrological modelling.Observations are needed to support more accurate parameterizations of biogenic CO 2 fluxes.We believe the bottom-up approach to model F C by SUEWS can be a promising tool in quantifying urban CO 2 emissions at the local scale, and therefore apply to capture the CO 2 emission hot spot, quantifying relative contribution of CO 2 sources, and assist to mitigate urban CO 2 emissions.

Appendix A: Workflow of LAI parameters optimization
In SUEWS, LAI influences the surface conductance, and subsequently Q E and F pho (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 matric 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 strong on "difficult functions" or larger dimensional search spaces (Hansen et al.).Taking Beijing as an example, the LAI parameters are optimized as the following: 1. City-level LAI derivation.The NASA Moderate Resolution Imaging Spectroradiometer (MODIS) product MOD15A2H with 8-day and 500-m resolution is first derived.The spatial average of LAI within the 6 th ring area of Beijing is extracted from the year 2008 to 2016, treated as the city-level of LAI.
2. Scaling city-level LAI to the tree-level.We note that the city-level LAI might be noticeably lower than the LAI at tree level due to the presence of non-vegetated surfaces in urban area.Nonetheless, the city-level LAI provides the signals 3. Spikes removal.Spikes occur occasionally due to instrument problems, uncertainties of retrieval algorithm and cloud contamination.Therefore, data points are identified as spikes if: where LAI 5 is the moving median of 5 consecutive data points, and β is constraint factor (here β = 2).The spikes are removed, marked as "sample discarded"; the other sample are kept, marked as "sample retained" and treated as "observed LAI".In the case shown here, 7 data points are discarded, accounting for 1.7% of the total sample.a obtained from a field survey over Beijing (Ma et al., 2019) b obtained by dividing maximum canopy conductance by G 1 (=3.5) The weighted average of g max is 7.0 mm s −1 for deciduous tree weighted by population ratio of each species (Table D1).
The average g max for grass is 3.7 mm s −1 .Note that when the species in the modelled area is known, we suggest the g max 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.
res vegetation and soil respiration.Positive values indicate sources of CO 2 and negative values sinks with respect to the atmosphere.

Figure 1 .
Figure 1.(a) The location of IAP tower and the land cover type within the 6 th Ring area of Beijing (Terra and Aqua combined Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover Type (MCD12Q1) Version 6) (Friedl and Sulla-Menashe, 2019), (b) a satellite image (Google Earth, image ©2022 Maxar Technologies) and (c) urban land use categories (EULUC-China) (Gong et al., 2020) covering the circle of 1 km radius around the IAP tower.
https://doi.org/10.5194/gmd-2022-305Preprint.Discussion started: 9 February 2023 c Author(s) 2023.CC BY 4.0 License.Traffic volume and fleet mix data for each link of the study domain for the year 2017 are obtained from Yang et al. (2019).

Figure 2 .
Figure 2. Annual average diurnal cycle of traffic rate (vehicle kilometers travelled, VKT) for weekday (VKT_WD), weekend (VKT_WE), and population density (POP) within the 1 km radius circle around the IAP tower.
Comparison of maximum conductances (gmax)  between the case base and case gs/gs_LAI.Parameters for case base followJärvi et al. (2011).More details can be found in Appendix D. gmax (mm s −1 ) case base case gs/gs_LAI to radius of modelled area SUEWS output will be evaluated against EC measured turbulent fluxes, but the challenge is that the source area of the observations are different to the exact modelling domain.To consider the impact of the chosen modelling domain on model evaluation, we designed three additional model runs where different radius circular areas around the IAP tower are considered.The default run is with the 1 km radius circle, but we also run SUEWS within 500 m, 750 m, and 1500 m circular areas, corresponding to https://doi.org/10.5194/gmd-2022-305Preprint.Discussion started: 9 February 2023 c Author(s) 2023.CC BY 4.0 License.withhigh R 2 (0.94) and low RMSE (0.4 m 2 m −2 ) (Appendix A).In the control case base, modelled LAI starts to rapidly increase from day of year (DOY) 70 and plateaus at DOY 105, which is too early when compared to remotely sensed LAI (MODIS LAI).Optimized LAI starts to grow at the same time but much slower and peaks around DOY 200.In autumn, LAI modelled by case base drops rapidly at DOY 310, while optimized LAI starts to decline steadily soon after its peak.LAI model with optimized parameters is better at capturing the behaviours of growth, peak growing season, and senescence for the year 2016 than in the control case base.

Figure 3 .
Figure 3. Normalized LAI in 2016, where the value 0 (1) represents the minimum (maximum) of LAI."Base" denotes the modelled LAI for the case base, "MODIS" the remotely sensed LAI within the 6 th Ring area of Beijing, and "optimized" the modelled LAI using the parameters derived with CMA-ES (Appendix A).

Figure 4 .
Figure 4. Input or modelled vs. observed hourly radiation fluxes, including (a-d) incoming solar radiation (K down ), (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 only K down is input and the rest are model output.

Figure 5 .
Figure 5. Average diurnal cycle of input or modelled and observed hourly radiation fluxes by season, including (a-d) incoming solar radiation (K down ), (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.Shading area denotes the standard deviation.Note that only K down is input and the rest are model output.

Figure 6 .
Figure 6.Annual and seasonal mean diurnal cycles of observed and modelled (a-e) latent heat flux (QE) and (f-j) sensible heat flux (QH ) for the four model runs (case base, gs, LAI, and gs_LAI) in the year 2016.

Figure 7 .
Figure 7. Model performance statistics R 2 , RMSE and MBE of (a-c) latent heat flux (QE) and (d-f) sensible heat flux (QH ) for the four model runs (case base, gs, LAI, and gs_LAI) in the year 2016.
photosynthesis rate is 5.0 µmol m −2 s −1 around noon in summer, while soil and vegetation respiration constantly serves as a CO 2 source with a rate lower than 3.6 µmol m −2 s −1 .Each of the model performance statistics of F C is of a similar magnitude among cases, indicating the F C modelling is far less sensitive to the choice of g max and LAI-related parameters than Q E and

Figure 8 .
Figure 8. Annual and seasonal average diurnal cycles of observed and modelled CO2 flux (F C ) for the four model runs (case base, gs, LAI, and gs_LAI) in the year 2016.

Figure 9 .
Figure 9. Seasonal average diurnal cycles of modelled CO2 flux (F C ) components by case gs_LAI in 2016.FcTraff denotes CO2 on-road traffic, FcBuilding building, FcMetab human metabolism, FcRespi vegetation and soil respiration, and FcPhoto vegetation photosynthesis.Positive values indicate sources of CO2 and negative values sinks with respect to the atmosphere.

Figure 10 .
Figure 10.Model performance statistics (a) R 2 , (b) RMSE and (c) MBE of CO2 flux (F C ) for the four cases (base, gs, LAI, and gs_LAI) in the year 2016.

Figure 11 .
Figure 11.(a) Land use fraction and mean daily traffic rate (VKT, same as T r in Table 1), and (b) accumulated CO2 flux (F C ) for modelling domain with the radius ranging from 500 m to 1500 m. "WD" denotes weekday, "WE" weekend; "Mod" denotes modelled F C , and "Obs" observed F C .Note that in (b) the lines for Mod 500m and Mod 750 m nearly overlap, and lines for Mod 1000m and Mod 1500m nearly overlap.
https://doi.org/10.5194/gmd-2022-305Preprint.Discussion started: 9 February 2023 c Author(s) 2023.CC BY 4.0 License.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 12 .
Figure 12.The contribution to total CO2 emissions by each component at the annual scale for the modelled area with a radius ranging from 500 m to 1500 m.Note that for photosynthesis the percentages denote the offset against total CO2 emissions.
https://doi.org/10.5194/gmd-2022-305Preprint.Discussion started: 9 February 2023 c Author(s) 2023.CC BY 4.0 License.ofvegetation phenology (e.g., leaf-out, peak growing season, leaf-fall).In order to give a more realistic estimate of LAI at the tree-level, the LAI city-level needs to be scaled.Therefore, LAI max and LAI min need to be given manually, preferably based on observational studies over local species.An observational study in Beijing shows that the maximum of LAI falls between 5 and 6 m 2 m −2(Wang et al., 2021).Here, LAI max = 6 and LAI min = 0.1 m 2 m −2 .
4. Interpolation.The observed LAI are linearly interpolated between values to obtain a daily time series, marked as "interpolated LAI". 5. Parameters derivation using CMA-ES.The time series of interpolated LAI and T air are subjected to CMA-ES to optimize the parameters.Using the LAI model and the parameters derived, LAI is calculated and marked as "predicted LAI".LAI model incorporated with the optimized parameters successfully reproduces the seasonal dynamics of observed LAI, capturing the timing of growth, peak growing season, and senescence (Fig. A1 a).The model performance is remarkably good, with high R 2 (0.94) and low RMSE (0.4 m 2 m −2 ) (Fig. A1 b).

Figure A1 .
Figure A1.(a) Time series of leaf area index (LAI) sample retained, LAI sample discarded, interpolated LAI and predicted LAI, and (b) comparison of the predicted against observed LAI.

Table 1 .
Ward et al. (2013) to anthropogenic heat and CO2 emissions.Minimum CO 2 release per capita 120 µmol CO 2 s −1 cap −1Ward et al. (2013) heat Fraction of fossil fuels used for heating 0.81 Cui et al. (2019); Zhang et al. (2019); MHURD (2018) F r nonheat Fraction of fossil fuels used for energy 0.2 this study E CO2perJ EF for fuels used in building energy use 0.1688 µmol CO 2 J −1 Cui et al. (2019) T base_HC Base temperature for heating degree day 12 • C this study Same as the case base but with g max values collected from observational studies over vegetation species in Beijing (Appendix D).The site-specific g max are in general lower than the values used by case base (Table3).4. Case gs_LAI: Same as the case base, but with both LAI and g max modified as described in case LAI and case gs.
Havu et al. (2022a)19)tions are: (1) parameters in the environmental response functions (Eq.4 -7) of surface conductance which follow those byHavu et al. (2022a)(Table2), because the product of response functions calculated followingKokkonen et al. (2019)is too low (95 th percentile = 0.19) to obtain realistic estimate of photosynthetic rate; (2) biogenic parameters for Eq.14-15 and soil properties which were updated followingHavu et al. (2022a)(Table2).In addition, F pho,max for grass/lawn (5.5 µmol m −2 s −1 ) is obtained from EC observations of CO 2 flux made over urban lawn in Helsinki in summer 2021 by fitting the conductance parameters (F pho,max,grass , G 2 − G 6 ) to the observations following Järvi et al.(2019) (Appendix C).2.Case LAI: Same as the case base but the parameters for Eq.1 describing the annual behaviour of LAI are optimized using remotely-sensed LAI and CMA-ES (Appendix A).The new optimized LAI parameters are compared with case base in Table4.3. Case gs: https://doi.org/10.5194/gmd-2022-305Preprint.Discussion started: 9 February 2023 c Author(s) 2023.CC BY 4.0 License.

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

Table 4 .
Järvi et al. (2011)area index (LAI) parameters between case base and case LAI/gs_LAI.All vegetation types (evergreen tree, deciduous tree and grass) use the same LAI parameters within one case.Case base values are as inJärvi et al. (2011).

Table 5 .
SUEWS model performance statistics for radiation fluxes.The abbreviations are the same as Fig. 4. Note that K down is an input and others are output of SUEWS.