Surface energy balance fluxes in a suburban of Beijing: energy partitioning variability

Abstract. Measurements of radiative and turbulent heat fluxes for 16 months in suburban Miyun with a mix of buildings and agriculture allows the changing role of these fluxes to be assessed. Daytime turbulent latent heat fluxes (QE) are largest in summer and smaller in winter, consistent with the net all wave radiation (Q*). Whereas, the daytime sensible heat flux (QH) is greatest in spring but smallest in summer, rather than winter as commonly observed in suburban areas. The results have larger seasonal variability in energy partitioning compared to previous suburban studies. Daytime energy partitioning is between: 0.15–0.57 for QH /Q* (mean summer = 0.16, winter = 0.46); 0.06–0.56 for QE / Q* (mean summer = 0.52, winter = 0.10), and 0.26–7.40 for QH  ⁄ QE (mean summer = 0.32; winter = 4.60). Compared to the literature for suburban aras, these are amongst the lowest and highest values. Results indicate that precipitation, irrigation, crop/vegetation growth activity and land use/cover all play critical roles in the energy partitioning. These results will help to enhance our understanding of surface–atmosphere energy exchanges over cities, and are critical to improving and evaluating urban canopy models needed to support integrated urban services, that include urban planning to mitigate the adverse effects of urban climate change.



Introduction
In the atmospheric boundary layer, turbulent flows are fundamental to mass and energy transport, so are 25 regarded as basic part of comprehensive observation studies of the urban boundary layer, such as the UBL/CLU-ESCOMPTE project in Marseille, France (Mestayer et al., 2005), the BUBBLE project in Basel, Switzerland (Rotach et al., 2005), and the SURF experiment in Beijing, China (Liang et al., 2018).
In this study, daily, monthly, and seasonal variations in surface energy fluxes for a suburban site (Miyun) in the Beijing are analyzed using 16-months of observations. We focus on the role of site 60 characteristics and environmental factors impacts on energy partitioning.

Site description
Miyun, located about 80 km northeast of the center of Beijing city (Tian'anmen Square) (Fig. 1a), has an area of 50 km 2 and population of ~500,000 (in 2019) (Beijing Miyun Statistical Yearbook, 2020). The
Within a 1 km radius of the Miyun instrumented tower, the surface is 69.6% impervious (buildings, roads, and parking lots) and 30.3% vegetation (wheat/maize rotation farmland, orchard, vegetable plots,

90
Both the eddy covariance (EC) system and radiometer are mounted 36 m above ground level (agl) on a 38-m triangular lattice tower, with the EC system pointing into the prevailing wind direction (east) and the radiometer south to avoid shadows. The EC system's three-dimensional sonic anemometerthermometer (CSAT3, Campbell Scientific Inc., USA) measures vertical, horizontal, and lateral wind velocity and virtual temperature; and the open-path infrared gas analyzer (LI-7500, LI-COR, Inc., USA) 95 measures water vapor and carbon dioxide molar densities. The 10 Hz data are logged on a CR3000 datalogger (Campbell Scientific Inc, USA). The 30 min turbulent sensible and latent heat fluxes are obtained using the EddyPro Advanced (v6.1.0 beta, LI-COR) software with standard correction procedures (Moncrieff et al., 1997) applied to ensure data quality (e.g., de-spiking raw data, tilt correction, time-lag compensation, double coordinate rotation, spectral corrections, and Webb et al. (1980) (Table 1). Data are excluded during rain and, 2 h after rain, or if low friction velocity (< 0.1 m s-1) occurred (Table 1). Wind direction data are corrected for changes in magnetic declination as EC sensor install based on magnetic north.
The radiation fluxes, measured with a CNR4 radiometer (Kipp & Zonen, Netherlands), are 1-min 105 samples by the CR3000 data-logger, from which the 30-min means are calculated. The incoming and outgoing longwave ( ↓ and ↑ ) and shortwave ( ↓ and ↑ ) radiation and net all-wave radiation ( * ) data are restricted to physically reasonable thresholds, with nocturnal shortwave radiation forced to 0 W m -2 (Michel et al., 2008). Additional observations analyzed are from an automatic weather station located 30 m from the EC tower, variables include: air temperature (HMP155, Vaisala, Finland) measured at 1.5 m agl, 115 precipitation from a tipping bucket rain gauge (SL3-1, Shanghai Meteorological Instrument Factory, China) mounted at 0.7 m agl (April to October), and from a weighing bucket rain gauge (DSC1, Aerospace Newsky Technology, China) mounted at 1.5 m agl (November to March); wind speed and direction (ZQZ-TF, Aerospace Newsky Technology, China) measured at 10 m agl, and soil temperature (ZQZ-TW, Aerospace Newsky Technology, China) measured at 0, 0.05, 0.10, 0.20, 0.40, 0.80, 1.6 and 120 3.2 m below ground. These data are sampled at 1 s, except for wind speed and direction (0.25 s), and averaged to 1 min (WUSH-BH, Aerospace Newsky Technology, China). Hourly means (or accumulated totals for precipitation) are used here. Quality control checks include: plausible range, internal consistency, temporal and spatial consistency, and other standard China Meteorological Administration network checks (Ren et al., 2015).  The study period precipitation differed significantly from Normal conditions in the region (Fig. 2b).
The 2012 was extremely wet, especially on a monthly basis with November (77.4 mm) and December (10.6 mm) about 580% and 400% greater than the Normal (13.4 mm, November; 2.7 mm, December), respectively. In contrast, most of 2013 had below average rainfall, except for June which was wetter 140 (186.7 mm, 223% of the Normal). Notably, dry spells occurred in May, November, and December 2013.
Only 0.7 mm (1.6% of Normal) rainfall fell in May and there was no rain in November and December 2013 at all.
Easterly winds (30-120º ) prevailed in MY every season in the night and for the whole day, with a frequency of 50% and 38% in spring, 54% and 43% in summer, 69% and 54% in autumn and 72% and 145 62% in winter . During the daytime, wind also comes from southwest direction (180-270º ) as a result of mountain-valley breeze, despite southwest wind differed in existence hours among seasons . In addition, the strongest winds (wind speed >5 m s -1 ) mainly came from the northeast (30-60º ), with a higher relative frequency in spring and winter ( Fig. 2c-  Footprint during unstable and neutral conditions are analyzed together by wind direction (5º sectors) 180 but split into day ( ↓ > 5 W m -2 ) and night. The median 90% source area extends to 582-677 m (by direction) during the daytime and 596-908 m at night (Fig. 1c), which corresponds to a daily median of 559-676 m. Within these source areas, the land use and land cover vary by sector from being more highly vegetated (30-150º ) to more built-up (210-360°) (Fig. 1d).

Data availability and analysis of fluxes
Missing data are not gap-filled. For analyses, that need 24 h from the mean fluxes are calculated before derived ratios (e.g. Albedo, Bowen ratio). Daytime and daily mean fluxes of net all-wave radiation, 205 sensible heat flux and latent heat flux are estimated based on monthly mean diurnal patterns.

Surface radiation budget
At the MY site, all radiation fluxes vary seasonally (Fig. 3). Daytime maxima of incoming short-wave radiation ( ↓ ) range from about 500 W m -2 in winter to 1000 W m -2 in summer ( Incoming longwave radiation ( ↓ ) is primarily influenced by near-surface air temperature and water vapor content (Flerchinger et al., 2009;Kotthaus and Grimmond, 2014a). Thus, the larger ↓ values are observed during warm and humid times of the year (i.e. summer -June-August) ( Fig. 3c) and smaller values observed in colder winter. The daily maxima of ↓ varies between ~320 W m -2 (winter) and ~470 230 W m -2 (summer). The monthly mean ↓ vary between 420 W m -2 (July) and 212 W m -2 (January 2013).
The latter was the coldest month of the observation period (Sect. 2.3).
Outgoing longwave radiation ( ↑ ) depends on the surface temperature and emissivity. The former is highly influenced by the total amount of incoming radiative energy. Thus, the larger ↑ is observed in summer ( Fig. 3d), consistent with the larger ↓ and ↓ . And vice versa, the lowest ↑ , along with smaller 235 ↓ and ↓ , measured in winter. The highest monthly mean ↑ of 467 W m -2 occurred in July and the smallest (275 W m -2 ) again in January 2013. Surface materials vary little over short periods, beyond the small impacts presence of snow and wet surfaces, but over multi-year periods external (e.g. building) materials will change having some influence. The latter is not relevant to this study.
As the net all-wave radiation ( * ) is net balance of the four radiation budget components, it is 240 affected by sun elevation angle, sky conditions, and surface characteristics. The variation of * is similar to that of the ↓ , with daytime maxima of * being greater in summer and less in winter. Daily mean * values vary between -20 W m -2 (winter) and 222 W m -2 (summer), with large scatter seen from May to https://doi.org/10.5194/acp-2021-1022 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License. September (Fig. 3e). Hence, the largest monthly mean * in summer, with the maximum value (126 W m -2 ) in July 2013. The smaller monthly mean * in June is attributed to its smaller ↓ caused by frequent 245 rainfall. Notably, the 4th June's extremely low ↓ results in a daily mean * than is < 0 W m -2 . With longer winter nights the monthly mean * is small or even negative. Notably, the minimum value during the observation period is in December 2012 when snow is present (i.e. increasing ↑ ) resulting in a monthly mean of only -0.5 W m -2 (Fig. 3e).  The median vary between 5 and 39 W m -2 diurnally across all of the months (Fig. 4). The fluxes are larger in summer and winter associated with cooling and heating needs, respectively. The monthly mean vary between 17 and 24 W m -2 (daily totals = 1.46-2.10 MJ d -1 m -2 ). These suburban values are larger than estimated in other suburban sites, such as 5 W m -2 in Basel (Christen and Vogt, 2004), 10-12 W m -260 2 in Montreal during winter (Bergeron and Strachan, 2012), and 6-10 W m -2 in Swindon (Ward et al., 2013). However, these suburban Beijing values are much less than the Beijing city-centre (~130 W m -2 , Wang et al., 2020). Hence, the MY values appear to be reasonable. The storage heat fluxes are determined using two methods (Sect. 2.6). Generally, storage heat flux is expected to have a net gain in summer and net loss in winter, giving almost zero net heat gain/release over the annual period (Grimmond et al., 1991). From this point of view, , ℎ has a credible annual variation in daily total values ( Fig. 5a) but may be a little low in winter. Although the , ℎ 270 underestimates more frequently during easterlies, given it is the prevailing wind direction there is no clear relation between underestimated , ℎ and wind direction or time of the day.
In autumn and winter (November -February) the differences in mean and median are negligible (< 0.5 W m -2 ) (Fig. 5). However, similar to residential Swindon results (Ward et al., 2013) , ℎ is smaller than , during the day but larger (more negative) at night (Fig. 5b). It is partly attributed to 275 our , ℎ coefficients not including all components of heat storage flux, such as biomass heat storage flux (Meyers and Hollinger, 2004;Oliphant et al., 2004). the relation between and * ) and 3 term (displacement of the and * peaks), whereas 2 is smaller (hysteresis; Grimmond et al., 1991;Oke and Meyn, 2009). These result in Swindon having a smaller diurnal range of , ℎ than at MY (Fig. S3).
As MY site has larger wintertime emissions, which direct impacts  )). This gives a lower limit, with , ℎ smaller than , − 2 in winter (Fig. S5). This indicates the winter , ℎ underestimates are related to OHM coefficients for the cropland wind directions, suggesting the wheat coefficients are poor (Fig. B2d). This is because smaller fluxes are obtained from using the soil heat flux ( ) data. These relations with * do not account for the crop biomass. Long-term datasets covering seasonal changes and 290 surface conditions, especially winter, can provide more OHM coefficients for use (Anandakumar 1999, Ward et al. 2013. , is used in the following analyses, but it is clearly biased to when data are available monthly daily total values are all positive throughout the year, even in winter (Fig. 5a).
In spring and summer, there is little differences in daytime values of , ℎ and , except for May 2013 (Fig. 5b). The May high , ℎ is the result of both fewer rain events and more radiative 295 input (Sect. 2.3 and 3.1; Fig. 2b, 3a). Under such conditions, more energy goes into heating the surface.
There is a sharp drop in , ℎ at night, mainly associated with a sudden decrease in soil surface temperature (0 m) (Fig. S4) that leads to the corresponding drop in * that impacts , ℎ (Appendix B). The rapid decrease in soil surface temperature is evident each month (around 18:00-19:00) with the transition from day to night (Fig. S4), and also with the valley wind to mountain wind (Fig. 1a, 2g-j).

300
The differences in soil temperature at night are less on clear day than on rainy day.   (Fig. 6a, b).

315
Sensible heat flux is greatest in spring (March-May) and smallest in June-September 2013 (Fig.   6c, d attributed to the extensive irrigated cropland in the prevailing wind direction (Fig. 1c, d; 2d, h). This enhances available energy to support evaporation and leads to smaller values (Dou et al., 2019).
Nocturnal is negative throughout the year, from a mean of -10.5 W m -2 in December (2013)

350
To facilitate comparison of surface energy balance (SEB) partitioning between sites a series of average daily, and daytime ( ↓ > 5 W m -2 ) flux ratios are analyzed. These include fluxes normalized by net allwave radiation ( * ) and the Bowen ratio (β= / ) (Fig. 7, 8; Table 2). The ratios are calculated at the 30 min time period allowing the distributions to be analyzed ( Fig. 7) but also using monthly mean fluxes.
During the observation period, daytime is 15-57% of * (Fig. 7). Unlike the absolute values 355 (Sect. 3.3), * ⁄ is largest in winter and smallest in summer. Whereas * ⁄ has the opposite shape with annual peak in July (Fig. 7). It varies between 6-56% during our study period.
At MY, from November to April (Fig. 7) dominates the * ratio, whereas from June to September dominates the * ratio. For May and October, they have a similar ratio with * (Fig. 7).  Oke et al., 1999). However, the * ⁄ ratios (0.08-0.18) are comparable to Tokyo ( * ⁄ = 0.07-0.11; Moriwaki and Kanda, 2004) and slightly larger than that of IAP Beijing ( * ⁄ = 0.07; Miao et al., 2012). This suggests the higher winter MY * ⁄ ratio could be attributed to its relatively lower heat absorbed and stored by buildings due to smaller (18.9 %). The is a small proportion of in the * energy partitioning.

Factors influencing energy balance fluxes
At MY, the impacts of precipitation on energy partitioning are obvious. As monthly rainfall from with vegetation cover of 2-31%) but still less than the values observed in central London (Kotthaus and   415 Grimmond 2012) or Marseille ).
Irrigation supplements precipitation and plays an important role in energy partitioning (Grimmond and Oke, 1986;Grimmond et al., 1996;Kokkonen et al., 2018). With only 3 days (26th, 27th and 28th May of 0.1, 0.4 and 0.2 mm, respectively) having rain in May 2013 (Fig. 2b), soil moisture decreases dramatically in May unless replenished by irrigation (Fig. 9). However, the May clearly decreased 420 because of the cropland irrigation (Fig. 10). The monthly mean of 0.50, clearly indicates the importance of this additional source of water.
Crop growth and phenology are also important influences. Transpiration during crop growth releases a large amount of water into the atmosphere observed as (e.g. Dou et al., 2019). The soil moisture in the irrigated cropland is higher than in the natural rainfed areas (Fig. 9). The gravimetric 425 samples dates are critical relative to the timing of both irrigation and rainfall (e.g. 8 July 2013 cf. few days later) as the soil moisture values can inverse, as cropland utilizes stored soil water via transpiration and evaporation.
In addition to these factors, land use/cover also plays a critical role in the energy budget. At MY, normalizing by the incoming radiation ( ↓ = ↓ + ↓ ) the ↓ ⁄ and ( ↓ ⁄ ) are always more (less) 430 when the wind is from directions with more buildings (cf. cropland-dominated) irrespective of season ( Fig. 11). The temporal variations in after rainfall differ between these two direction types (Fig. 12).
Water can drain relatively quickly after rain in the building-dominated direction because of the high proportion of impervious surface. Whereas, water can infiltrate into the soil (i.e. stored) so can remain near 1 for longer periods after rainfall (Fig. 12).    Table 1

450
In this analysis of surface energy flux measurements for a suburb (MY) of Beijing over 16 months we gain a better understanding of surface-atmosphere dynamics.
All five components of the radiation balance ( ↓ , ↑ , ↓ , ↑ , * ) have unimodal annual patterns, with higher values in summer (lower in winter) with the maximum monthly means in July (minima December / January).

455
At MY, daytime sensible heat flux is greatest in spring. However, it is smallest in summer rather than winter, unlike previous suburban studies. This is a because of the prevailing wind direction has extensive cropland, irrigation, and frequent rainfall in summer. All of these factors are thought to play a role. Nocturnal is negative throughout the year and smaller in winter, so the minimum daily total still occurs in winter.

460
The latent heat flux is positive throughout the day. The daytime maximum is greatest in summer July (lowest in winter December) because of the influence of rainfall, irrigation, plant growth activity and available energy. Monthly median diurnal maxima of vary from 10-248 W m -2 during measurement period. These are close to both the smallest and largest values reported for suburban areas.
The maximum monthly daily total is in summer (July 2013) and minimum winter (December 2013). traffic , , and human metabolism , (Allen et al., 2011;Lindberg et al. 2013). In this study, we 485 determine appropriate temperature response coefficient for the MY area.

A1 Temperature response coefficient
In LQF, the mathematical expression of , is as follows: where is the population density (capita ha -1 ) and , is daily building energy consumption 490 (kWh day -1 per capita). The latter, from daily electricity consumption data, accounts for 16% of the total energy consumption in MY (Table A1). Several key parameters to determine are shown in Fig. A1. Here is the air temperature when 500 the energy consumption is the lowest. If the air temperature is higher (lower) than , more energy will be consumed due to cooling (heating). ( ℎ ) is the building energy consumption thermal response slope for cooling when air temperature above (below) , also known as a cooling (heating) coefficient.
is the minimum energy consumption. and are threshold values for air temperature. When air temperature is beyond this range, energy consumption has reached saturation and no longer increases 505 with air temperature changes.
A new set of temperature response parameters applicable to this region is obtained by using daily air temperature and electricity consumption data for the Miyun district in 2012-2013 (Table A1). The resulting parameters are given in Table A2. The energy consumption response to temperature changes is a "V"-shape curve in Miyun, rather than a "U"-shape as seen in Shanghai city (Ao et al., 2018). The 510 cooling coefficient in MY is lower than that in Shanghai city (0.04 º C -1 ), which is attributed to relatively short period air conditioning in MY given the cooler summer nights caused by topography and continental monsoon climate. However, in winter with colder air temperatures and a longer heating period the heating coefficient ℎ for MY is higher than for Shanghai ( ℎ =0.01º C -1 ) (Ao et al., 2018).

A3 Population density
The where is time and, i is the area fraction covered by th land cover type. Given the differences in land cover around MY flux tower, the plan area fractions by 30° wind sector footprint area are used to 545 calculate ∆ , ℎ by direction (Fig. B1).
The coefficients used (Table B1) include values obtained from 30 min measured data of net radiation * and soil heat flux in Weishan Farmland Experimental Station (116° 09′ E, 36° 39′ N). Weishan and MY have similar climate background, the same wheat/maize rotation pattern, and consistent growth cycle (Lei and Yang, 2010;Lei et al., 2018).

550
The Weishan station has a 10 m tower with a four-component radiometer (CNR1, Kipp & Zonen, Netherlands) installed at 3.5 m agl and two soil heat flux plates (HFP01SC, Hukseflux, Netherlands) installed at a depth of -0.03 m. The instrument data are recorded at 10-min intervals using a CR10X data logger (Campbell Scientific, USA). The soil heat flux data are the average of two measurement sites on the east and west of the tower.

555
For application convenience, the cropland coefficients are determined by season. A random, but arbitrary 70% of data are selected in each season to fit the coefficients. The remaining 30% is used to https://doi.org/10.5194/acp-2021-1022 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License.
evaluate the simulation results. The results are good, but poorest in winter (Fig. B2). These values are used for the farmland (Table B1).
Building coefficients also vary with season (Table B1). Here we use data from Nanjing city (Wang 560 et al., 2008), with an average of summer and winter values used in spring and autumn. Other coefficients are used year round. For the Road/Impervious areas an average of literature values is used (Table B1).   The observation data used in the study are available from the corresponding author with permission (Email: jxdou@ium.cn).

Author contribution
JD conducted the eddy covariance measurement, carried out the data analyses, and wrote the draft. SG supervised the scientific interpretation of the results and polished the writing. SM, BH, HL and ML 575 provided key data sets and contributed to the anthropogenic heat flux estimation, and OHM coefficients fitting.