Spatial and temporal patterns of surface-atmosphere energy exchange in a dense urban environment using scintillometry

Spatially integrated measurements of the surface energy balance (SEB) are needed in urban areas to evaluate urban climate models and satellite observations. Scintillometers allow observations of sensible heat ﬂux ( Q H ) over much larger areas than techniques such as eddy covariance (EC), however methods are needed to partition between remaining unmeasured SEB terms. This is the ﬁrst study to use observed spatial and temporal patterns of Q H from a scintillometer network to constrain estimates of remaining SEB terms in a dense, heterogeneous urban environment. Results show that Q H dominates the surface energy balance in central London throughout the year, with expected diurnal courses and seasonal trends in Q H magnitude related to solar radiation input. Measurements also reveal a clear anthropogenic component of Q H with winter (summer) weekday Q H values 11.7% (5.1%) higher than weekends. Spatially, Q H magnitude is correlated with vegetation and building landcover fraction in the measurement source areas. Spatial analysis provides additional evidence of anthropogenic inﬂuence with highest weekday/weekend ratios (1.55) from the City of London. Spatial differences are used to estimate horizontal advection and a novel method to estimate monthly latent heat ﬂux is developed based on observed landcover and wet–dry surface variations in normalized Q H . Annual anthropogenic heat emissions are estimated to be 46.3 W m − 2 using an energy balance residual approach. The methods presented here have potential to signiﬁcantly enhance understanding of urban areas, particularly in areas with tall buildings where there are few observational data.


Introduction
Urban areas are growing rapidly and contain the majority of the world's population.Urbanization physically modifies thermal and radiative properties of the surface and drastically alters surface-atmosphere energy exchanges in cities. Intense human activity in urban areas also creates anthropogenic heat emissions which can be an important energy source to the surface energy balance.The partitioning of surface-atmosphere energy exchange is a fundamental control on urban climate processes and is required for many applications such as predicting urban air temperatures, simulating urban boundary-layer dynamics, and modelling air pollutant transport and dispersion.To improve understanding of urban energy balance processes, there is a need for spatially integrated measurements to test urban climate models and compare with satellite observations.

Urban surface energy balance
The surface energy balance for an urban area can be expressed as: where Q * is net all-wave radiation, Q F is anthropogenic heat input, Q H is turbulent sensible heat flux, Q E is turbulent latent heat flux, Q S is net heat storage in the three-dimensional urban surface, and Q A is net heat advection.Here, we use sign convention where positive Q * and Q S represent downward fluxes (i.e.surface energy gain), and positive Q E and Q H are upward fluxes (i.e.surface energy loss).Positive Q F is an energy input and positive (negative) Q A represents a net energy gain (loss) to the system.Units are W m −2 .
The Q A term is typically assumed to be negligible if there is sufficiently homogenous upwind fetch.However, in coastal cities Q A can be important due to land-sea thermal circulations (i.e.sea breezes) with values exceeding 100 W m −2 and decreasing with the distance inshore (Pigeon et al., 2007).
The storage term, Q S , is important on hourly to weekly timescales, but tends to zero over annual time-scales (e.g.Offerle et al., 2005).On daily scales, buildings and pavement in the city typically warm throughout the day and act as a storage reservoir of thermal energy.This energy is either lost via long-wave radiative cooling or acts to warm the air (where it is measured primarily as Q H ).This heat release can persist overnight and drive unstable nighttime boundary layers (Grimmond and Oke, 1999b).Buildings can also act as a storage buffer for Q F and increase time between heat emissions inside a building and subsequent release to the atmosphere (Bohnenstengel et al., 2014).
Anthropogenic heat emissions (Q F ) are a unique feature of urban environments and result from energy consumption (fuel combustion, air conditioning/heating, electrical appliances/lighting) and human metabolism (Sailor, 2011).Q F varies on diurnal, weekly, and seasonal time-scales according to human activity patterns (e.g.commuter traffic cycles) and energy demands.Q F is difficult to measure directly on neighborhood to city scales and there are several approaches that have been used to quantify its magnitude.The residual approach observes all other terms of the urban energy balance and since Q F is measured as part of the other terms (including Q * via long-wave radiation losses), the difference between the sum of all terms and Q * is taken to be Q F (Offerle et al., 2005).Other approaches include inventory-based methods that rely on fuel and electricity consumption statistics or modelling approaches such as up-scaling building energy model outputs (Heiple and Sailor, 2008).
Another notable feature of the urban energy balance is modified energy partitioning towards higher Q H and lower Q E (i.e. higher Bowen ratio), compared to more natural surroundings (Oke, 1988).This is largely due to relative lack of surface moisture and vegetation transpiration in cities and has implications for urban boundary-layer heights and thermal turbulence production.Q H is also influenced by Q F .Iamarino et al. (2012) found the majority of Q F in London is expressed as Q H (80.7%), with smaller portions going into Q E (7.3%) and wastewater (12.0%).Modelling work by Bohnenstengel et al. (2014) indicates that about two-thirds of Q F goes into Q H with the remaining one-third serving to warm the surface and increase outgoing long-wave radiation losses.

Scintillometer measurements in urban areas
The focus of this article is on observations of Q H given its importance to the surface energy balance and role in urban boundary-layer dynamics.Observations of Q H using scintillometers have advantages over other measurement techniques, such as eddy covariance (EC), because measurements are integrated over a path length of several kilometres instead of the sensor being located at a single point.As a result, scintillometer measurement source areas can be much larger than those of EC measurements.
Most studies have compared scintillometers with EC results and found generally good agreement between the two techniques, especially during daytime unstable conditions.Several sites have reported scintillometer fluxes which tend to slightly overestimate EC fluxes, in particular for low values (i.e.overnight) (e.g.Lagouarde et al., 2006;Zieli ński et al., 2013).During night-time and day-night transition periods (stable, near neutral), there is generally more scatter in EC-scintillometer comparisons because measurements may be above the surface layer (SL) and Monin-Obukhov similarity theory (MOST) assumptions begin to break down (Ward et al., 2015a).
In Gongju, Lee et al. (2015) compared scintillometerderived sensible heat and momentum fluxes to outputs from a high-resolution (200 m) numerical model configured for their heterogeneous urban setting.Results showed good agreement between the model and scintillometer methods with estimated uncertainty in scintillometer Q H of 20-30%.
In general, urban scintillometer studies report a smoother temporal signal than EC, attributed to larger spatial integration of observations.Magnitudes of Q H have all been within expected ranges based on values reported in the urban energy balance literature from EC observations.In Ł ódź and Swindon, monthly variations in Q H (lower magnitudes in winter, higher in summer) are observed based on variations in short-wave and net radiation.In Swindon, Q H values are linked to vegetation landcover fraction (negative correlation) through source area modelling.This study also reported reductions in Q H following rain events by analyzing high-resolution (10 min average) temporal signals.
The majority of these sites have been relatively densely built (e.g.built fraction = 61% in Tokyo, 58% in Marseille).Exceptions are the suburban residential areas of Swindon (vegetated fraction = 53%) and Helsinki.The latter two sites are notable both for the distribution of surface types integrated by observations and the distribution of observed atmospheric stability.There is a greater frequency of near-neutral and stable conditions relative to the more densely built sites (attributable to reduced storage heat releases and anthropogenic heating overnight) which affects Q H magnitude and choices related to flux calculations (section 2.3).
When deriving Q H from scintillometer measurements, stability can be used to determine whether use of specific free convection assumptions is appropriate.Comparison between the 'free convection' and 'mixed convection' solutions to similarity theory functions to calculate Q H indicate good agreement during unstable conditions (Lagouarde et al., 2006;Zieli ński et al., 2013).Overall, the free convection solution tends to underestimate relative to the full solution during near-neutral conditions due to neglect of mechanical turbulence (section 2.3).However, an advantage of the free convection solution is that it relies on fewer meteorological variables and associated uncertainties (Lagouarde et al., 2006).
Together, these studies have demonstrated the feasability of using scintillometers to measure Q H in cities, however there remain challenges in interpreting results over heterogeneous urban surfaces.In particular, there is a need to develop methods to constrain estimates of partitioning between the remaining unmeasured energy balance terms in Eq. ( 1).There is also need to explicitly account for spatial variations during data analysis and to develop methods to exploit this variability when possible.The broad objective of this research is to investigate the spatial and temporal patterns of urban Q H using a unique scintillometer network configuration in central London.Findings are then used to develop new methods to partition between energy balance terms on monthly and annual time-scales.

Site description
A network of three scintillometers has been established to measure area-averaged Q H in central London, UK (Figure 1).2).Source area calculations are described in section 2.5.P2 is used as exemplary for all paths.
provides complete annual coverage with continuous winter data (December-February).
Average monthly mean air temperature during the 12 month study period observed at London Heathrow Airport was 11.63 • C (average daily maximum = 16.05 • C, average daily minimum = 8.24 • C), average annual precipitation was 561.2 mm, and average hours of sunshine totalled 1485.2 h (17.0%) (UKMO, 2016).The mean temperature for the study period is in the 79th percentile and total precipitation is in the 35th percentile of all continuous 12 month periods since measurements began in 1948.
Total sunshine hours are in the 36th percentile of all 12 month periods since measurements began in 1958.
Physically, buildings and impervious surfaces (concrete and asphalt roads, pavement, etc.) dominate the central London study area and combined they account for nearly 80% of plan area landcover fraction (section 2.5, Figure 2).Mean building height is around 15 m (section 2.5, Figure 2) and building stock is diverse with 37% of Inner London buildings constructed before 1900 and 11.5% built after 2000 (GLA, 2016).Buildings house nearly 30 000 m 2 of commercial, retail, and office floorspace in Inner London, which is 23.5% of all office floorspace in England (GLA, 2016).There are also several large vegetated parks in the region (Hyde Park/Kensington Gardens ∼4 km 2 ; Regent's Park ∼2 km 2 ; Green Park/ St. James's Park/Buckingham Palace Gardens ∼1 km 2 ; Figure 1).The study area is also intensely modified by daily and weekly cycles of anthropogenic activity.In 2014, the population of Inner London was 3.4 million, with 233 292 residents in the City of Westminster and 8072 in the City of London (GLA, 2016).These two boroughs are highlighted because they make up a significant portion of the scintillometer source area (along with Islington) and are generally representative of central London in terms of anthropogenic activity (Figure 1).Though the residential population of these two areas is less than 10% of total Inner London population, 1.16 million people pass into central London in the morning rush hours between 0600 and 0900 on an average weekday (Monday-Friday) (Björkegren et al., 2015).In Westminster and City of London combined, there were over 1 billion vehicle kilometres travelled by all vehicles in 2014 (GLA, 2016) despite the presence of a congestion charge zone (in operation from Monday to Friday 0700-1800 local time) (Figure 1) which has reduced vehicle kilometres travelled by 23% between 2000and 2012(TfL, 2016)).Tourism is also a significant driver of anthropogenic activity in central London.In 2014, there were over 17.4 million international visits to London, an average of nearly 48 000 visits per day (GLA, 2016).
This concentration of employment, commuter, and tourist activity alters patterns of energy consumption and anthropogenic heat release compared to more residential areas.For example, in Westminster and the City of London combined, annual energy intensity (defined as total natural gas and electricity consumption divided by total floor space) of commercial areas is 697 MWh m −2 .For residential areas, energy intensity is only 1.7 MWh m −2 (GLA, 2016).Total commercial energy consumption (natural gas + electricity) is nearly six times greater than residential energy consumption in these boroughs.

Instrument configuration
The scintillometer network is composed of two LASs (LASMkII, Kipp & Zonen, Delft, The Netherlands) and one boundary-layer scintillometer (BLS) (BLS900, Scintec, Rottenburg, Germany).The LAS instruments both use a single infrared beam (wavelength 850 nm) and the BLS operates with two parallel infrared beams separated by 0.164 m (wavelength 880 nm).The instruments have a receiver aperture diameter of 0.15 m.
The sensors are oriented so that their beams form a triangle (Figure 1, Table 1).Path 1 (P1) is between Islington Michael Cliffe House (IMU) and the BT Tower (BTT) (BLS sensor); Path 2 (P2) is between BTT and Barbican Cromwell Tower (BCT) (LAS sensor); and Path 3 (P3) is between BCT and IMU (LAS sensor).Mounting heights are 91 m asl (above sea level) at IMU, 180 m asl at BTT, and 145 m asl at BCT.These heights are used to ensure that observations are above the blending height of surface roughness elements (2-3 times mean building height in urban areas) and are within the atmospheric surface layer (usually defined as the lowest 10% of the boundarylayer height, z i ) so that application of surface-layer similarity theory to determine Q H is valid.Ceilometer measurements in central London show boundary-layer heights regularly grow to 2000 m during mid-day in summer and usually exceed 500 m even in winter, which indicates scintillometer measurements are within the surface layer during these times.An additional advantage of elevated measurements is that signal saturation is reduced because the temperature structure parameter decreases with height.
The effective beam height (z fb ) of each scintillometer path length (L p ) is determined using the average height of the beam above the ground surface weighted by the scintillometer path weight function at 1 m horizontal resolution along the beam (Hartogensis et al., 2003).Surface elevations are from a LiDAR digital elevation dataset flown during 2008 (Lindberg and Grimmond, 2011).The L p calculation also takes into account vertical displacement between transmitter and receiver.The effective measurement height (z f ) is calculated as z f = z fb − z d where z d is wind direction-dependent source area weighted displacement length (section 2.5).
Prior to the scintillometer deployment in the current configuration, instruments were run simultaneously along a single path length at a different location in London (L p = 1290 m, z f = 51 m) (Mustchin et al., 2013).Results from this intercomparison show that 1 min structure parameter (C 2 n ) values from the two LAS sensors are in good agreement during a representative summer 24 h period (P2 LAS versus P3 LAS: slope = 0.91, R 2 = 0.97).During a second 24 h period, a single LAS sensor (P3 LAS) was compared against the BLS and 1 min C 2 n values are also in good agreement (P3 LAS versus P1 BLS: slope = 0.89, R 2 = 0.98).
Auxiliary meteorological data from the London Urban Meteorological Observatory (LUMO, 2016)  velocities, virtual acoustic air temperature, and water vapour and carbon dioxide mixing ratios were recorded at 10 Hz and flux processing was performed in 30 min intervals.Detailed descriptions of the EC site and data processing are given in Kotthaus andGrimmond (2012, 2014a,b).Additional details and analysis of F CO2 observations are given in Ward et al. (2015b) and Björkegren et al. (2015).
Incoming (K↓) and outgoing (K↑) solar radiation, and incoming (L↓) and outgoing (L↑) long-wave radiation are also measured at KCL at a height of 60.9 m asl (CNR4, Kipp & Zonen) and allow for calculation of For precipitation, a spatial radar rainfall dataset from the UK Met Office available for 1 km × 1 km grid cells at 5 min resolution is used (UKMO, 2003).Nine grid cells (3 km ×3 km total area) centred on the scintillometer configuration were selected as representative of rainfall within the scintillometer source areas.The radar datasets are verified against rain-gauge data by the Met Office to minimize error and comparisons of accumulated totals indicate the radar product tends to overestimate precipitation compared to gauges, depending on the radar range (Harrison et al., 2000).For this analysis, the radar product is used only to indicate whether it was raining or not; the total amount is not used.

LAS
This section describes the calculations and assumptions to determine Q H from raw scintillometer measurements of beam intensity (I).The general steps are: (i) calculate path-averaged structure parameter (C 2 n ) using I according to wave propagation theory, (ii) calculate the temperature structure parameter (C 2 T ) using C 2 n , and (iii) derive Q H from C 2 T using surface-layer similarity theory.
At the receiver, I varies in time due to pertubations in the scintillometer beam resulting from air density variations across the path length.The density variations are spatially integrated resulting in path-averaged C 2 n , calculated for a 1 min averaging period (Wang et al., 1978): where the instrument-specific coefficient α = 1.12,A is the variance of the natural logarithm of I fluctuations (σ 2 lnI ) during the 1 min averaging period, D is the LAS aperture diameter, and L p is the scintillometer path length between the transmitter and receiver.
The C 2 n values can be defined in terms of individual structure parameters for temperature (C 2 T ), humidity (C 2 q ), and the covariant term (C Tq ) (Hill et al., 1980): where A T and A q are structure parameter functions (Andreas, 1988), T is mean air temperature, and q is mean humidity.This formulation neglects pressure fluctuations as they occur over time-scales longer than the 1 min averaging time and on spatial scales greater than L p .Further simplifications can be introduced because the infrared beam wavelength is primarily sensitive to temperature fluctuations, so moisture-related fluctuations are neglected.An additional assumption is that of a generally dry surface based on EC energy balance measurements at KCL with a mean Bowen ratio of 5.1 during the 12 month study period.This assumption makes the correction to account for the influence of water vapour fluctuations negligible.Therefore, when surface conditions are dry, Q H is then calculated from C 2 T MOST using similarity functions for unstable conditions: (5) and stable conditions: Although there are several similarity functions in the literature, those of Wyngaard et al. (1971) are selected here given their frequent use in scintillometer literature and that they meet MOST theoretical requirements of a neutral limit and −2/3 scaling (Kooijmans and Hartogensis, 2016).
Two sets of c 1−3 coefficients are compared in this article to assess Q H uncertainty: those of Andreas (1988) (An; c 1 = 4.9, c 2 = 6.1, and c 3 = 2.2) and De Bruin et al. (1993) (Db; c 1 = 4.9, c 2 = 9, and c 3 = 0).These coefficients are chosen because they are routinely applied in urban scintillometer studies (e.g.Lagouarde et al., 2006;Zieli ński et al., 2013;Ward et al., 2014) so provide a means to assess uncertainty in our results.However, several sets of coefficients are found in the literature (e.g.Li et al., 2012).Recently, Kooijmans and Hartogensis (2016) explored MOST parameters using data from 11 scintillometer sites and covering a range of stability conditions, and recommend c 1 = 5.6 and c 2 = 6.5 during unstable conditions.Unfortunately, this very useful study contains no urban datasets.Research into MOST parameters in urban areas is an active area of research and urbanspecific coefficients have been developed for use in the roughness sub-layer using small-aperture scintillometers by Kanda et al. (2002) and Roth et al. (2006).Ward et al. (2015a) examined six sets of MOST coefficients in Swindon, UK and opted for those of Andreas (1988).It is beyond the scope of the current study to explore these important research questions concerning MOST coefficients in urban areas in depth; we merely aim to quantify uncertainty at this site related to choice of MOST coefficients and maintain compatibility with other urban sites by using the An and Db coefficients.
The stability parameter, ζ , is defined as (z fb − z d )/L, where z d varies by 10 • wind direction sectors (section 2.5) and L is the Obukhov length: where u * is friction velocity, g is gravitational acceleration, and k is von Kármán's constant (0.4).The temperature scaling variable, T * , is calculated as: Friction velocity is estimated using a stability-adjusted logarithmic profile and dynamic roughness length (z 0 varies based on 10 • wind direction sectors, section 2.5).Measured wind speed at the scintillometer transceiver is also adjusted to z fb using the logarithmic wind profile.The wind profile is then solved iteratively with Eqs ( 5)-( 8) and the sensible heat flux is calculated as: where ρ is air density and c p is specific heat capacity of air.During unstable free convection conditions when ζ → −∞, a simplified MOST solution for Q H is (De Bruin et al., 1995): where b = 0.57 (b = c 2 ).This solution is advantageous because it depends on fewer meteorological inputs (and associated uncertainties) compared to the full solution (MX) and explicit choice of Q H sign is not required during day-night stability transitions because Q H is assumed always positive.However, the free convection solution (FC) takes into account only buoyancygenerated turbulence and ignores mechanical production and thus underestimates Q H when free convection conditions are not met (De Bruin et al., 1995).

BLS
The BLS instrument includes pre-processing corrections (Scintec, 2016) that are not applied to LAS data.During conditions of high turbulence, strong scattering of the scintillometer beam results in a saturation effect where further increases in turbulence are not captured by I measurements.The BLS uses a numerical look-up table based on Clifford et al. (1974) to correct for saturation effects.Mean saturation correction for the BLS is < 1% and data during periods when the saturation correction is greater than 5% are excluded from further analysis.
A second correction is applied to account for absorption fluctuations.These fluctuations are caused by processes other than turbulent mixing, such as fog or dust, which occur on spatial length-scales on the order of the instrument path length.The BLS is able to correct for these fluctuations because of the dual-beam configuration and path separation of the instrument (Scintec, 2016).The annual mean absorption fluctuation correction is −3.1%.
Both corrections are applied during calculation of the σ 2 lnI signal.The remainder of flux processing then follows Eqs ( 2)-( 10), where in Eq. ( 2) α = 4.629 and where σ I is the variance of I and I is the mean of I (Scintec, 2016).

Quality control
Several quality control filters are implemented, both on raw measurements and processed fluxes (Figure 3(a)).First, to remove spurious datapoints in I, logical thresholds (min = 0, max ≈ 3× daily maximum) are applied and data are withheld during instrument maintenance periods.Next, data are withheld during rain events (measured by a tipping-bucket rain gauge at BCT) because the sensor cannot measure I properly when there is precipitation in the scintillometer path.Similarly, fog occludes the signal and datapoints where I is less than 0.5 times the daily maximum value are used to indicate fog conditions and are withheld (Ward et al., 2014).This step also removes rainy periods that the rain gauges fail to register due to spatial and temporal resolution limitations of the mechanical gauges.Saturation is a recognized issue that affects scintillometer measurements during certain atmospheric conditions and instrument configurations (Kohsiek et al., 2006).This occurs when the level of scintillation surpasses a structure parameter (C 2 n ) threshold (S th ) and further increases in scintillation become de-coupled from observed changes in I.The S th is a function of L p , aperture diameter D, and instrument wavelength (λ) (Kleissl et al., 2010): When the S th value is exceeded, LAS C 2 n values are withheld from further processing.This step is not applied to the BLS data because there is onboard saturation correction software (section 2.3.2).
Observations during stable atmospheric conditions are also deemed unreliable because the scintillometer beams are likely above the atmospheric surface layer (i.e.lowest 10% of the boundary-layer height z i ) and assumptions to calculate fluxes using similarity theory begin to fail.Specifically, (i) the observation height is higher than 0.1×z i and measurements are de-coupled from the surface flux and (ii) vertical divergence of Q H within the surface layer can be high from entrainment during periods of boundary-layer growth (Braam et al., 2012).To filter these periods, data when air temperatures measured at two height levels (z m1 = 145 m asl, z m2 = 21 m asl) at BCT indicate stable conditions (T zm1 − T zm2 > 0 K) are withheld.A second step uses a night-time (K↓<1 W m −2 ) maximum Q H threshold of 2.5 times the observed daily minimum.Values above this threshold are presumed unrealistic at night and indicate atmospheric conditions are not conducive for measurements.This step is designed to filter stable conditions at larger spatial scales which may not be registered by the urban micro-and local-scale temperature profile measurements.

Source area calculation
Flux source areas are calculated using the Kormann and Meixner (2001) source area model calculated on a 14.4 × 14.4 km domain at 4 × 4 m resolution.A statistical look-up table approach is used based on observed frequency of atmospheric conditions (instead of calculating a source area for every time step) to minimize computational expense.In all, 644 source areas were calculated at the centre point of each scintillometer beam using various combinations of binned input parameters.
The binned input parameters are wind direction, L, u * , and standard deviation of horizontal wind velocity (σ v ).Source areas are calculated for every 10 • using combinations of L, u * , and σ v .Sensitivity tests of these inputs on the size of the 50% cumulative source area weighting extent were used to determine appropriate bin size.Three L inputs are used to represent stable (L = 10 m), unstable (L = −10 m), and near-neutral (L = 1000 m) stability regimes.Source area dimensions are not overly sensitive to L explicitly, but are sensitive to u * and σ v , which vary with stability.The u * is binned in increments of 0.25 m s −1 from 0.25 to 1.50 m s −1 .The σ v and u * are dependent variables, so for each u * value a σ v value is calculated from a linear relation between u * and σ v determined from sonic anemometer observations at KCL during the study period (slope = 1.7, y intercept = 0.17).Additional σ v values ±0.3 m s −1 from the central value determined from linear regression are also used (in total three σ v values for each u * value).
The z 0 and z d inputs to the source area model are determined using an iterative procedure.First, an initial estimate is performed based on average building height (z b ) within the source area model domain.Source areas are modelled using these inputs and new z 0 and z d values are determined based on source areaweighted z b , where z 0 = 0.1z b and z d = 0.7z b (Grimmond and Oke, 1999a).This iteration is performed three times before source area-weighted mean z 0 and z d values for 10 • wind direction bins are used as input for Q H calculations (section 2.3).
So far, this procedure produces source areas calculated at the L p centre-point at the effective beam height.To determine the source area influencing the entire L p , the centre-point source area is replicated along the pathlength of the scintillometer beam horizontally every 50 m and multiplied by the scintillometer pathweighting at that point.This procedure introduces uncertainty because the height of the beam and z 0 and z d of the upwind area influencing each point along the beam vary.The uncertainty is mitigated by the path-weighting function which is most heavily weighted towards the beam centre-point.
The annual cumulative source area is determined as the weighted mean of the 644 individual source areas (Figure 1).Source area weights are determined by the observed frequency at the KCL sonic anemometer of wind direction and atmospheric conditions (cumulative frequency is 90%).This procedure assumes identical conditions at the sonic anemometer (60.9 m asl) and across the scintillometer L p .This is justified theoretically because measurements from EC and the scintillometers are usually situated in the constant-flux layer above the roughness sublayer.
This same procedure is used to calculate the annual cumulative source area for the KCL EC station (Figure 1) and differs from methods used in previous studies for that site (e.g.Kotthaus and Grimmond, 2014b).
Source area-weighted landcover fraction (f ) is determined by multiplying annual source area weights (φ) for each 4 × 4 m grid cell (x, y) by the plan-area landcover of the modelling domain.The resulting matrix is then summed for each landcover classification (i) (road, building, soil, water, deciduous vegetation, coniferous vegetation, grass): On average, less than 100% of the source area falls within the modelling domain (Table 2).For the area outside the domain, landcover proportions are assumed to be the average of the portion within the domain.
To test sensitivity of landcover fraction to the source area calculation method, the path-integrated source area for each path was calculated using two methods: (i) the method described previously using the centre-point calculated source area replicated at 50 m increments along the path length, and (ii) individual source areas using point-specific z fb , z 0 , and z d at 50 m increments along the path length.
Identical source area model meteorological inputs representative of convective conditions (u * = 1.25 m s −1 , σ v = 2.3 m s −1 ) from the primary wind direction (230-240 • ) are used for each method for all paths.These conditions are chosen because unstable conditions have generally smaller source areas than stable conditions and thus landcover fractions are expected to be more sensitive to source area uncertainties.
Next, landcover fractions (f i ) for each path are calculated using Eq. ( 13) for both methods.For all landcover classes, differences in f i between methods are < ±1% for all paths.For P1 and P2, method (i) produces a slightly larger source area (1.75 versus 3.38 km 2 and 1.42 versus 2.70 km 2 cumulative 50% weighted area).For P3, method (ii) produces a larger path-integrated source area than that of method (i) (1.27 versus 1.56 km 2 cumulative 50% weighted area).

Q H intercomparisons and sensitivity
This section describes a set of three comparisons to assess the scintillometer measurement and data processing techniques: (i) scintillometer measurements are compared directly against simultaneous EC measurements, (ii) Q H calculation techniques are compared, and (iii) flux calculations are compared for each path.Finally, the sensitivity of flux processing to input variables is tested.For these comparisons, summer (JJA) data are used because this season has highest data recovery rates.

Comparing eddy covariance and scintillometers
First, scintillometer measurements are directly compared to EC measurements from KCL to check if scintillometer values are reasonable.Perfect agreement is not expected because scintillometers and EC sample different source areas (Figure 1).For this analysis, 1 min scintillometer measurements are averaged to 30 min intervals to match temporal resolution of EC measurements.Then Q H data are sorted by K↓ into day (K↓>1 W m −2 ) and night periods (K↓<1 W m −2 ).
During daytime, overall agreement between 30 min scintillometer and EC Q H is moderate (Figure 4(a)).At lower values (Q H < 200 W m −2 ), the scintillometer tends to measure higher Q H than EC, though measurements converge at larger values (Q H 400 W m −2 ).There is significant scatter between measurement techniques (mean R 2 = 0.50), though this is to be expected because of differences in measurement source areas.
Overnight, agreement between 30 min scintillometry and EC is not as good as during daytime (Figure 4(a)).This is not unexpected given the theoretical limitations of both measurement techniques during night-time operation (e.g.vertical flux divergence in urban canopy layer as the atmosphere approaches near-neutral or stable conditions) and is consistent with other comparisons in the urban scintillometry literature (e.g. Ward et al., 2014).For all sites and solutions, the scintillometers measure higher Q H relative to EC when EC measurements are low (Q H < 100 W m −2 ), and lower Q H relative to EC when EC measurements are high (Q H > 100 W m −2 ) (mean slope = 0.13).
Nocturnal Q H in central London is maintained by Q F and Q S releases and Q H > 100 W m −2 is not unusual.
A portion of the EC-scintillometer difference may be explained by the proportion of water present in the measurement source area.The EC tower is located near the River Thames and the long-term EC turbulent flux source area contains 12.1% water surface cover, compared to 1.9-3.7% in the scintillometer source areas (Figures 1 and 2).The river's surface temperature is affected by tidal flow patterns, but is expected to be cooler than the overlying air during summer daytime, on average, and thus act as a local heat sink (negative Q H ).There is also expected to be greater energy partitioning toward Q E from the river (i.e.lower β).Previous EC studies at KCL have found the river's influence to be greatest when winds are aligned parallel with the river upwind of the EC sensors and that channelling flow above the recessed riverbed (approximately by 5-10 m) is likely to result in local-scale advection (Kotthaus and Grimmond, 2014b).
During summer (JJA), ensemble hourly median Q H from scintillometers (FC solution) are 26.1% higher than EC throughout the 24 h diurnal period, though agreement overnight between ensemble means for all sensors is good (Figure 5).To investigate if this difference can be explained by the river's presence, observed Q H for all sensors is scaled to the nonwater surface of the respective flux source area.This assumes Q H = 0 W m −2 over the river on average, though overnight the river is likely warmer than the overlying air and thus would support positive Q H and local instability.After this scaling, the difference between EC and scintillometers is reduced to 18.5% (relative to scintillometers) throughout the day on average.This correction may be slightly overestimated (if Q H over the river is positive overnight), yet still is not enough to explain a large portion of the EC-scintillometer difference.
Another likely reason for EC underestimation of Q H relative to scintillometers is due to lack of energy balance closure for EC systems.EC typically underestimates turbulent heat fluxes on the order of 20% in forest and other natural ecosystems due to several causes, including undersampling of highfrequency turbulent eddies, insufficient flux averaging periods, and advection (Wilson et al., 2002).Large-eddy simulation studies have suggested secondary circulations caused by heterogeneous landscape interactions with the flow are a primary cause of energy balance non-closure (Kanda et al., 2004;Inagaki et al., 2006).Scintillometers appear to improve energy balance closure due to their ability to spatially average over larger areas and measure energy transport by secondary circulations that EC measurements at a single point cannot resolve (Foken et al., 2010;Liu et al., 2011).
As a first-order estimate to investigate if EB closure can explain the EC-scintillometer differences (after scaling to take into account the river's influence), a typical 20% energy balance closure residual is assumed for the EC system (e.g.Wilson et al., 2002) and the scintillometer is assumed to measure Q H perfectly (for now, this sets aside scintillometer errors).If the residual is added proportionally to the EC-measured Q H , scintillometer-EC Q H difference is reduced to 10.4%, on average.

Comparing scintillometer Q H calculations
Next, the MX and FC flux calculation solutions are compared for each sensor path (Figure 4(b)).During daytime, agreement between the free convection solution (FC) and mixed convection (MX) with de Bruin coefficients (Db) is very good at all sites.This is to be expected as the b coefficient in the FC solution is based on the Db coefficients (b = c 2 ).The MX equation with Andreas (An) coefficients tends to underestimate FC and MX-Db by 9 to 15%.At night, FC and MX-Db are in very good agreement at P2, whereas MX-An underestimates both FC and MX-Db.At P3 and P1, both MX-An and FC are in agreement and tend to underestimate MX-Db.Note that stable periods have been removed during quality control (section 2.4).
Good comparison results between MX and FC solution are achieved because free convection conditions are often met.This is primarily due to the frequent occurrence of small L values and to a lesser extent the high effective measurement height of the scintillometers.To assess the impact of measurement height on distribution of the stability function ζ , ζ is determined using the KCL measurement height of 50.3 m agl and L calculated from measurements at the KCL EC tower.From this, 10.8% of available 30 min periods are classified as near-neutral or stable (ζ > −0.01).Next, ζ is calculated using the effective beam height of the scintillometers (98.7, 111.6, 140.7 m agl, Table 1) with identical L from KCL.Using these heights, the proportion of near-neutral ζ is reduced to 8.0% for P1, 7.6% for P2, and 8.3% for P3.
In summary, FC and MX-Db flux calculation solutions agree very well and MX-An tends to give lowest fluxes during the day (−9 to −15%).In general, at night there is more scatter in all comparisons.

Comparing scintillometer paths
The third comparison is between scintillometer paths, for each flux calculation solution (Figure 4(c)).During daytime, agreement between sites is generally very good for all solutions, particularly between P1 and P3.However, P2 has a constant offset of 28.5 W m −2 relative to P1 and 42.2 W m −2 relative to P3.At night, the same pattern is observed with P2 higher relative to P1 and P3.Given the nearly constant nature of this offset (Figure 4(c)), a systematic bias resulting from uncertainties in z f and L p seems the most likely explanation (Figure 6); however, differences between sensors should be somewhat mitigated because all path pairs share an endpoint in the sensor configuration (section 2.4, Figure 1).
There are several additional factors that contribute to uncertainty in comparisons between sensors.Flux divergence below measurement height is likely a factor when comparing measurements at different heights.Source area differences between paths are another source of uncertainty and it is possible that the P2 source area includes a higher proportion of the central business district with higher Q H and Q F .Observed  path differences are also likely caused by differences between instruments, though this source of uncertainty is expected to be minor (section 2.2).

Q H sensivitity to input variables
Next, the sensitivity of Q H calculations to meteorological and surface morphometry input variables is assessed (Figure 6).This analysis was conducted for P3 during a midday period representative of summertime midweek, convective conditions (n = 240).Input variables tested are: z 0 , z d , z f , wind speed, air temperature, pressure, and L p .Data for this period were processed using Db coefficients with a single input variable perturbed and compared to a baseline using unperturbed input variables.It should be noted that z 0 , z d , z f are not independent because in this study building height determines z 0 and z d , which in turn affects z f .
For both solutions (MX and FC), Q H is most sensitive to z f with a change of ±5 m resulting in change of fluxes of ±6.3-6.5% (Figure 6).The MX solution is also sensitive to wind speed (±3.5% with ±1 m s −1 change), however the FC solution does not use wind speed as an input.This is consistent with results of a similar sensitivity analysis by Hartogensis et al. (2003).
An additional source of uncertainty is the choice of empirical calculation coefficients.Choosing the An instead of Db coefficients in the mixed convection solution alters Q H by −13% and use of Kohsiek (Ko) (Kohsiek, 1982) coefficients instead of Db changes results by −5% for the FC solution.This is also illustrated in the linear comparison between Mx-Db and Mx-An, where regression slopes range from 0.85 to 0.90 (Figure 4   Although there is uncertainty with both solutions, results using the FC assumption will be analyzed in the remainder of the article.This is not intended to be applicable during stable or neutral periods.Fortunately, independent temperature profile and sonic anemometer measurements indicate there are very few periods (<4%) with stable conditions and these periods have been removed from analysis during quality control (section 2.4, Figure 3).Sonic anemometer measurements at KCL indicate near-neutral conditions are also infrequent (|ζ | < 0.01 ≈ 4%).Additionally, the FC solution requires fewer input variables (and associated uncertainties) because it does not rely on the logarithmic wind profile.
In summary, scintillometer measurements tend to measure higher Q H than EC measurements at this site, especially during daytime.This is consistent with observations at other urban sites in Marseille (Lagouarde et al., 2006) andŁ ódź (Zieliński et al., 2013).If all techniques and assumptions for MOST scintillometry and EC are met perfectly, scintillometer-EC differences would be due entirely to source area differences.In reality, differences are due to a combination of factors, including: (i) uncertainties regarding similarity theory coefficients and meteorological variables used for scintillometer Q H calculations, (ii) uncertainties determining instrument measurement height and corresponding flux divergence below measurement height, (iii) source area differences between sensors (in particular the influence of the River Thames), and (iv) relative lack of EC energy balance closure.Scintillometer sites agree well with each other during the day, though P2 has a constant positive offset compared to P1 (28.5 W m −2 ) and P3 (42.2 W m −2 ), likely from uncertainties in determining z f and flux divergence below measurement height.

Q H temporal patterns
Diurnal patterns of Q H follow an expected pattern based on incoming solar radiation and net radiation (Figure 5).During summer, median hourly overnight (2300-0500 h) (all times in the remainder of the analysis are UTC) values of Q H are positive on average (50 W m −2 for all paths), indicating an unstable atmosphere.Net radiation at KCL is negative for these hours (mean = −50 W m −2 ), implying Q S and Q F emissions maintain Q H .
After sunrise, Q H rises steadily and peaks at 1300 h, 1 h after K↓ and Q * peak.Throughout the afternoon and early evening (1400-2000 h), a portion of storage heat releases and anthropogenic heat emissions are manifested as Q H , resulting in an asymmetric decline of Q H relative to the morning rise.This hysteresis pattern between Q * and Q H is consistent with other observations in the literature (e.g. Ward et al., 2014).
The mean midday (1000-1400 h) ratio of Q H :Q * for all scintillometers during JJA is 0.70.This Q H :Q * ratio is higher than observed in other European cities such as Basel (0.66 during SON and 0.46-0.49during JJA 1100-1500 h) and Ł ódź (0.34-0.56 during Aug-Sep daytime).These urban areas had slightly more vegetation landcover (0.16-0.31 in Basel and 0.24-0.78 in Ł ódź) than the present site (0.16-0.19), which acts to decrease the Q H :Q * ratio.A source of uncertainty in the ratio calculated in this study is the spatial discrepancy between the radiometer source area and scintillometer source area.Additionally, use of the Db MOST coefficients tends to overestimate Q H relative to An and Ko coefficients (section 3.1).
Scintillometer-measured Q H magnitude varies seasonally relative to the annual mean for all paths through the 24 h day (Figure 7).The annual mean is determined as the equally weighted mean of the ensemble diurnal course from all four seasons.During spring (MAM) and summer (JJA) overnight hours (2300-0500 h), Q H is 0-25% higher for all paths than the annual average.During midday (1000-1800 h), spring and summer measurements are 25-70% times higher than the annual average.Autumn (SON) and winter (DJF) Q H are below the annual mean throughout the 24 h diurnal course.Overnight, Q H is 0-50% below the annual mean, throughout the day Q H is 50% below average, and greatest departure (50-75%) below average occurs during 1700-1900 h.Greatest departures from the annual mean (both positive and negative) occur during afternoon and early evening  which suggests Q S release from the urban surface regulates Q H at these times of day.
There are also differences between Q H measured during weekdays (1800 h Sunday-1800 h Friday) and weekends (1800 h Friday-1800 h Sunday and public holidays, i.e. non-working days) (Figure 8).These cut-off times were selected because Friday and Saturday evenings are 'weekend' nights and Sunday evening is functionally a 'weekday' night for many people in London.For this analysis, individual weekdays and weekend days were selected with similar amounts of incoming short-wave radiation (daily cumulative totals within 25%) to control for Q H differences from solar energy input.Weekday-weekend differences are expressed as ratios for each sensor individually to avoid biases caused by systematic differences between sensors (section 3.1).
During summer, weekday Q H at 0900 is 18.0% higher than on weekends and at 1600 is 23.0%higher.This two-peaked curve resembles the diurnal pattern of weekday commuter vehicle traffic and associated heat emissions found in many cities (e.g.Sailor and Lu, 2004).In central London, vehicle traffic generally remains high throughout the day with no obvious secondary peak (Ward et al., 2015b), however scintillometer source areas are likely influenced by areas outside the city core with traffic patterns more representative of suburban regions (Figure 1).Early evening (2000-2200 h) weekday Q H is higher than weekends, but overnight (2300-0300 h) weekday Q H is on average less (−15.9%)than weekend Q H , likely due to shifts in anthropogenic activity and traffic patterns.Daily total weekday Q H is 5.1% greater than weekend daily totals during summer.
In winter, the diurnal course of the weekday:weekend Q H ratio differs from that in summer.Overnight values are similar to summer (−20% relative to weekends), but a large increase (44.8%) is observed relative to weekends at 1000 h.Weekday daytime values remain above weekend values until 2100 h.The likely explanation for weekday-weekend differences is Q F emissions from indoor space heating during the working week.Earlier sunset and later sunrise during winter also means more energy is consumed and heat released for indoor and outdoor lighting.For winter, daily total weekday Q H is 11.7% greater than at the weekend.

Q H spatial patterns
The distribution of scintillometer-measured Q H magnitude over central London is influenced by the distribution of wind directions and atmospheric stability via modification of turbulent flux source area orientation and shape.Over a homogoneous surface (i.e.length-scale of surface variations length-scale of measurements), changes in source area shape and orientation would be irrelevant because the surface composition of heat sources and sinks remain the same.However, in urban environments surface variations of landcover, heat sources, and sinks can be of the same order of magnitude as length-scales of measurements, so variable source areas must be taken into account.
For this site, mean source area weighted vegetation landcover fraction (f V ) is highest (25.1%) in the northwest quadrant (270-360 • ) and northeast quadrant (25.7%) (0-90 • ) for all sensor paths (Figure 2(a)).Lowest f V (11.3%) is found in the southeast quadrant (90-180 • ).Building landcover fraction (f B ) follows the inverse pattern with highest f B in the southeast and southwest (41.8 and 43.5%, respectively) and lowest in the northwest and northeast (32.0 and 31.3%).Water landcover (f W ) is a minor component (mean = 2.9%) and is largest when winds are from the south (up to 9.2% for 180-190 • ) and source areas include a portion of the River Thames.Impervious fraction (f I ) is generally constant in all wind directions (mean = 41.1%).
Mean building height (z b ) is 14.3 m and is highest in the southeast (17.4 m) and southwest (16.4 m) and lowest in notheast (11.1 m) and norhtwest (12.1 m) (Figure 2(a)).The z 0 and z d also follow the same spatial pattern as z b because they are directly scaled to z b in this study (z 0 = 0.1z b and z d = 0.7z b ; section 2.5).
To analyze spatial variability, Q H is binned by wind direction (20 • increments) and Q H magnitude (100 W m −2 bins) during summer (JJA) (Figure 9(a)).The most frequently observed wind direction at BCT when Q H data are available is between 220 • and 240 • (21.1%).Wind directions with highest proportion (30.9%) of Q H >300 W m −2 are from 140 • to 160 • .In contrast, wind directions with highest proportion (91.8%) of Q H < 300 W m −2 are from 0 • to 20 • .
To isolate spatial differences in Q H , both anthropogenic influences and synoptic conditions (e.g.rain, variations in K↓) must be controlled for.To do this, Q H is conditionally sampled during warm-month (May-September) working days (nonholiday Monday-Friday) when 200 <K↓< 400 W m −2 and it has been at least 15 h since the most recent rain event as measured by radar.(Here, a rain event is defined as precipitation in any of the nine 1 × 1 km 2 radar grid cells; section 2.2).The 15 h threshold is based on findings from EC measurements at KCL showing the Bowen ratio is reduced up to 12 h after a rain event (Kotthaus and Grimmond, 2014a).These controls are meant to ensure that comparison between wind directions are from times with a dry surface, similar solar radiation input, similar levels of anthropogenic activity, and that sufficient number of data points are retained for analysis (n > 10 for each wind direction sector).
The selected data are then binned in 20 • sectors (w) and binned median values ( Q H ) are normalized (Q H * ) as the ratio of individual bin medians ( Q Hw ) to the equally weighted mean of binned medians ( Q H ) from all wind direction sectors: This procedure expresses spatial differences as ratios for each sensor individually to avoid biases from systematic differences between sensors (section 3.1).
On average, highest weekday Q H * from all sensors is measured when winds are from the south (100-260 • ).Q H * is 14% higher on average than Q H when winds are from 100 • to 270 • and are 11% lower than Q H when winds are from the north (280 • to 80 • ).Highest mean Q H * from all sensors occurs when winds are from 180 • to 200 • (19.5% higher).Flux source areas with  weekend Q H ratio of hourly medians for all seasons (box plots), winter (solid line) and summer (dashed line).P1 is excluded during winter due to insufficient data because of instrument maintenance (Figure 3), and the dashed line is mean weekday:weekend ratio for JJA for all sensors.
southerly winds generally include high density areas of buildings, traffic, and anthropogenic activity.
When winds are from the northwest (280-360 • ), Q H * is 12% below average (Figure 9(b)), likely due to the increased presence of vegetation and moisture in Regent's Park (Figures 1 and 2(a)).When winds are from the northeast (0-80 • ), Q H * is 10% below Q H on average due to nearly equal amounts of f V and f B from 10 • to 50 • (Figure 2(a)).
An anthropogenic signal is seen in the wind direction analysis as well (Figure 9(c)).For this analysis, weekday (working day) observations are compared with weekend (non-working day) observations.All other environmental controls (K↓, dry surface, warm weather months) are identical to those in Figure 9(b).On average for all wind direction sectors and sensors, weekday fluxes are 22.6% higher than at the weekend.Mean weekend-weekday differences are smallest (−3%) in the  14)) in 20 • wind direction bins for all scintillometer sites, and (c) weekday:weekend Q H ratio of binned medians in 15 • bins.For (b)-(c), section 3.3 gives details on conditional sampling criteria to control for season, day of week, solar radiation input, and dry surface conditions.
northwest sector (260-360 • ).Weekend-weekday differences are highest when winds are from 80 • to 120 • ) (55% above the spatial average) and 160-220 • (51% above average).When winds are from 80 • to 180 • , source areas include the City of London, a predominantly commercial and financial centre with a high proportion of office buildings (Figure 1).During the working week, daytime population in this area grows to 553 103; nearly 70 times greater than the weekend/overnight residential population of 8072 (GLA, 2016).When source area-weighted landcover fractions are plotted against Q H * in 20 • bins (Figure 10) during dry surface conditions (identical controls as in Figures 9(b) and (c)), both vegetation and building landcover fraction are linearly correlated with Q H * .As the landcover fraction of vegetation increases (and buildings decrease), Q H * also decreases due to greater partitioning of net available energy toward evapotranspiration and reduced contributions from Q F and Q S .This relation is consistent with results in the urban energy balance literature from other cities (e.g.Grimmond and Oke, 2002;Christen and Vogt, 2004;Loridan and Grimmond, 2012).
It is then possible to analyze Q H * immediately following precipitation when the surface is wet (i.e.<3 h since the most recent rain event with precipitation in all nine 1 × 1 km 2 radar grid cells) (Figure 10).Q H * for wet periods is defined relative to Q H for dry periods, as in Eq. ( 14).All other controls remain identical to those used in Figures 9(b) and (c).The 3 h threshold is somewhat arbitrary, though findings from EC measurements at KCL show the Bowen ratio is reduced up to 12 h after a rain event (Kotthaus and Grimmond, 2014a).Although the surface is 'wet', the Bowen ratio measured by EC remains high (median β = 3.1 during JJA from 1000 to 1800 h and time since rain <3 h, compared to β = 5.3 during the same period when time since rain >15 h) so that scintillometer flux calculation assumptions that neglect humidity remain valid (section 2.3).The Bowen ratio remains high after rain (though reduced compared to dry surface state) due to anthropogenic heat releases, efficient chanelling and drainage of water from the surface before it can be evaporated, and often significant areas of wall that remain dry from rain shadow effects.
On average, Q H * magnitude is reduced by 27% for all wind direction sectors when the surface is wet (Figure 10).Somewhat surprisingly, vegetation (building) fraction and Q H * remain negatively (positively) correlated for both wet and dry periods when Q H * for wet surface conditions is normalized by the wet Q H .It was hypothesized that rain would act to reduce horizontal thermal variations between landcover types and that Q H would be more spatially homogenous during wet periods (i.e. the linear slope would be reduced).This effect was demonstrated in a modelling study by Ramamurthy and Bou-Zeid (2014).A possible explanation for the observed pattern is that the correlation between landcover fraction and Q H is maintained by the presence of vertical walls that remain dry and greater mechanically generated turbulence from areas with high building landcover.

Advection of vertical sensible heat flux
Horizontal gradients in Q H (Figure 9(b)) in the presence of mean flow signifies horizontal advection ( Q A ) of vertical Q H .A firstorder estimate of intra-urban Q A values at this site is determined using observed horizontal Q H differences ( Q H ) between wind direction sectors in Figure 9(b) along with mean mixed-layer wind speeds (U) and horizontal length-scale of scintillometer source areas ( x) for a time step (t): This approach is conceptually a single-layer box-model with two adjacent, well-mixed columns of equal mixed-layer height but differing vertical Q H .For this analysis, t is set to 1800 s (typical 30 min flux averaging period), U is set to a representative mixed-layer velocity of 5 m s −1 (based on observations at KCL and BCT), and x is set to 5000 m as a typical horizontal distance of the scintillometer source area 80% cumulative weighting area (Figure 1).The weekday Q H values are determined for each wind sector in Figure 9(b) as , where w is the centre-point ( • ) of the wind sector bin and Q Hw is the mean of all sensors.
In this scenario, daytime weekday values of Q A are up to ±134.4 W m −2 when winds are along a southwest-northwest axis ( Q H = ±74.7 W m −2 ).Here, Q S is positive when flow is from relatively high Q H to low Q H (from southwest to northwest) and negative when flow is from low to high (from northwest to southwest) (Figure 9(b)).Smallest Q A values (±20.8W m −2 ) occur when flow is along an east-west axis where Q H = ±11.6W m −2 .These values of Q A are in line with those calculated in Marseille using observations of horizontal air temperature gradient (Pigeon et al., 2007).
This analysis applies only to daytime advection of vertical Q H within the long-term flux source area, however, and does not consider larger-scale processes such as thermal (e.g.sea breeze, urban-rural breeze) or synoptic circulations.Similarly, energy advection associated with horizontal gradients in moisture is unaccounted for.

Latent heat flux
Latent heat flux in an urban setting is composed of evaporation (E), transpiration (T ), and anthropogenic contributions (Q FE ) from sources such as fossil fuel combustion, respiration, and releases from air conditioning: Here, we estimate monthly cumulative Q E using measured Q H and a modified Bowen ratio approach.Monthly Q E is modelled using precipitation data and relations between measured Q H and landcover developed in section 3.3.In this approach ('landcover method'), T , E, and Q FE are considered separately: where E D and E W represent bulk evaporation during times when the surface is dry and wet, F D and F W are dry and wet surface fractions at monthly time-scales, τ is a coefficient representing vegetation and soil evapotranspiration, and f V is vegetation landcover fraction.
Collectively, the bracketed terms in Eq. ( 17) are approximately 1/β, but this excludes the Q FE term.Although the E D and E W terms are conceptually similar to β terms for wet and dry conditions, they only account for evaporation, not transpiration.In this equation, Q H is the cumulative total of hourly ensemble median Q H for each month (Figure 11).Overall, Q FE is expected to be a minor contribution to Q E .Modelling work in London calculates this to be 7.3% of total Q F , or approximately 1-10 W m −2 on average in greater London (Iamarino et al., 2012).Since this term is minor compared to measured vertical fluxes and their uncertainties, we neglect Q FE in this analysis.
First, evapotranspiration from vegetated surfaces is calculated using τ derived from the slope of the linear relation (τ = 1.22) between vegetation fraction and normalized Q H during dry surface conditions (Figure 10).Using this slope and the source area weighted f V (9.6%) from the southwest primary wind direction quadrant (180-270 • ), T is set to 0.117 Q H .This makes the simplifiying assumptions that linear relations developed during warm-weather daytime periods apply also to daily energy exchange totals throughout the year.During leaf-off months (roughly October-April), this introduces uncertainty because grass and coniferous trees transpire at a reduced rate relative to summer, and deciduous trees are largely inactive primarily due to reductions in photosynthetically active radiation (Kramer and Boyer, 1995).Additional biological and environmental processes driving plant transpiration such as vapour pressure deficit, wind speed, vegetation type and health, and stomatal resistance are obviously also not included in this method.This estimate also includes an evaporation component from soil beneath aboveground vegetation.These uncertainties are somewhat mitigated by the overall small fraction of vegetation in the source area and thus presumably small contribution of transpiration to Q E (Table 2).
Next, the evaporation coefficient E W is calculated from observed differences in Q H when the surface is wet compared to dry (Figure 10).On average for all sensors and all wind directions, Q H is reduced by 27% when the surface is wet (Figure 10) and this procedure assumes this fraction of energy is partitioned entirely towards evaporation.Thus, E W is set to 0.27 and E D is set to 0, which assumes minimal evaporation from exposed water sources when the surface is dry.Again, this is conceptually similar to a Bowen ratio approach because we assume energy is partioned towards the evaporation component of Q E when the surface is wet.This method introduces uncertainty through use of a binary wet or dry surface, when in reality Q H changes continuously after a rainfall event as surface water evaporates.EC measurements at KCL indicate Q H is reduced for up to 12 h after precipitation (Kotthaus and Grimmond, 2014a).
The F D and F W fractions are estimated using radar precipitation data collected for nine 1 × 1 km 2 grid cells (section 2.2).For each month, F W is calculated as the number of days where the radar records precipitation in all nine cells divided by the total number of days in the month (NDays rain /NDays total ).All other times are classified as dry (F D = 1 − F W ).
For comparison, independent Q E measurements from the KCL EC station are also used and a third estimate of Q E is made using the Bowen ratio (Eq.( 16)) for each month measured by the KCL EC system (Figure 11).The EC Bowen ratio for each month is calculated as the 24 h sum of median EC Q H divided by the 24 h sum of median EC Q E .Q E for each month is then calculated using the monthly 24 h sum of hourly median scintillometer Q H and 1/β (Eq.( 16)).These methods are not expected to agree exactly, given the source area differences between EC and scintillometer measurements, however the comparison provides a useful reference.
Using the landcover method, monthly Q E ranges from 7.6 W m −2 (mean daily total 0.66 MJ m −2 day −1 ) in November to 27.2 W m −2 (mean daily total 2.35 MJ m −2 day −1 ) in July (Figure 11).Mean monthly β for the entire year using the landcover method is 5.2, with values ranging from 4.3 in December to 6.2 in October.The annual course of monthly Q E totals agrees well with independent EC measurements and annual cumulative Q E is 1.3% lower than the EC measurements (Figure 11).The monthly courses of the Bowen ratio and landcover method are also in good agreement, though annual cumulative β Q E overestimates the landcover method by 13.4%.
An advantage of the landcover method is that Q E is scaled to Q H and can be implemented with easily available precipitation data and landcover information.This technique would be useful for estimating energy balance partitioning in situations where Q H is already available, for example if Q H is calculated from remotely sensed surface temperatures.However, given the uncertainties involved, it is only recommended for low-vegetation landcover environments where Q E is a relatively minor component of the surface energy balance.

Residual, storage, and anthropogenic heat
Observed Q H and modelled Q E are used to constrain the amount of energy available for Q F and Q S using an energy balance residual approach (Figure 11).For each month, hourly ensemble medians of Q H measured by the P3 scintillometer, Q E modelled using the landcover method, and Q * measured at KCL are calculated.The P3 scintillometer data are used because it has the highest data recovery rate for all months (Figure 3).Hourly medians are summed across the 24 h diurnal cycle and residual flux (Q R ) is calculated as The Q R values also include the portion of Q F expressed as outgoing long-wave radiation loss due to warming of the 3D urban surface (which acts to reduce Q * ).This is expected to be approximately one-third of total Q F based on model results by Bohnenstengel et al. (2014).The method does not account for Q F that is transported out of the study area via wastewater, which has been modelled to be approximately 10% of total Q F in London (Iamarino et al., 2012).In general, there is uncertainty with any residual approach because all measurement and modelling errors are accumulated in Q R , therefore Q R is very sensitive to small errors in the other flux terms.Additional uncertainty arises because the radiative flux source area of the radiometer does not match the turbulent flux source area of the scintillometers, though the radiative source area is assumed to be representative of the larger turbulent flux source area.
Q R follows an annual cycle with highest Q R value during February (61.3 W m −2 , daily total 5.3 MJ m −2 day −1 ) and lowest during June (30.1W m −2 , daily total 2.6 MJ m −2 day −1 ).Monthly Q R is linearly correlated with monthly total heating degree days and Q R (Eq. ( 18)) and annual mean Q F .Q H is measured at P3, Q * is measured at KCL and inner-decile ranges are shown as shaded regions.Three values of Q E are shown: (i) 'Q E landcover' is modelled using the landcover method described in section 3.5, (ii) 'Q E EC' measured at KCL by EC, and (iii) 'Q E β' is estimated using the EC Bowen ratio and P3 Q H (section 3.5).The Q E shaded region is the range across all three values for each month.Q R is determined from Eq. ( 17) and the Q R shaded region represents the Q R range using the maximum and minimum of values of the shaded Q H , Q * , and Q E regions.Section 3.6 gives an explanation of annual mean Q F .
(HDD) (Figure 12).HDD are calculated using the difference between observed daily mean canopy layer air temperature within the scintillometer source area at 3 m agl near the P3 reciever location and a UK-specific reference temperature of 15.5 • C (Day, 2006).
Q R is also linearly correlated with net carbon dioxide emissions (F CO2 ) measured at KCL (Figure 12).Q F and F CO2 are closely linked because heat and CO 2 emissions are generated primarily from fossil fuel combustion and released by human metabolic activity.This suggests that Q R at this location is dominated by Q F which tends to be driven largely by natural gas combustion for indoor space heating (e.g.Iamarino et al., 2012).Uncertainty with this comparison is from heat sources that do not generate local CO 2 emissions (e.g.waste heat from electrical appliances and lighting), CO 2 emissions from soil and vegetation respiration, CO 2 sequestration by vegetation, and Q S .
There is a hysteresis evident in both the HDD Q R and F CO2 Q R relations likely due to Q S .The residuals between Q R and the linear models range from −8.8 to 8.9 W m −2 (−0.76 to 0.77 MJ m −2 day −1 ).These numbers are reasonable bounds of monthly Q S and consistent with results in central European cities.Based on observations in Basel (Christen and Vogt, 2004) and Ł ódź (Offerle et al., 2005), monthly Q S is expected to be within ±23.1 W m −2 (±2 MJ m −2 day −1 ) with net positive Q S (i.e.energy uptake by the surface so that Q R < Q F ) from January to June and net negative Q S (i.e.surface energy losses so that Q R > Q F ) from July to December.
If Q S is 0 W m −2 over the course of a full year (i.e.no net energy gain or loss), then annual cumulative Q R is equal to Q F (Figure 11).Annual cumulative Q R is determined as the sum of monthly Q R (Eq. ( 18)) multiplied by the number of days in each month.Distributed evenly throughout the year, average annual Q F is 4.0 MJ m −2 day −1 , or 46.3 W m −2 .This is comparable to Q F results from other climatically similar cities, including modelling studies in London.Iamarino et al. (2012) modelled mean annual Q F for the entire Greater London region as 10.9 W m −2 and for the City of London as around 140 W m −2 for 2005-2008.In Chicago and Philadelphia, city-scale inventories determined peak summer Q F values range from 30 to 60 W m −2 and from 70 to 75 W m −2 during winter (Sailor and Lu, 2004).In Ł ódź, Q F averaged 32 W m −2 from October to March (Offerle et al., 2005) and in Basel annual mean Q F was 20 W m −2 using an energy balance residual approach.

Summary and conclusions
This work presents a unique dataset from a network of three scintillometers in central London for a full year (December 2014-November 2015).It is the first study to use observed temporal and spatial patterns of Q H to estimate partitioning between remaining energy balance terms Q E , Q A , and annual Q F .The novel methods developed here offer significant new potential for extracting additional information about surface-atmosphere energy exchange processes from observations.
Temporally, Q H exhibits expected diurnal and seasonal cycles related to solar insolation.There is also a strong anthropogenic signal in measured Q H . Winter weekday cumulative Q H is 11.7% higher than weekends from building energy occupancy and in summer, weekday cumulative totals are 5.1% higher than weekends.A diurnal commuter traffic signal is detectable in both winter and summer.
Observed Q H also varies spatially according to wind direction when conditional sampling is used to control for anthropogenic and synoptic factors.Normalized Q H varies with vegetation and building landcover fraction in the scintillometer measurement source area.Weekday Q H is also found to be up to 55% higher than weekend Q H when wind directions are from the central business district (the City of London).Spatial differences are then used to estimate the magnitude of horizontal advection of Q H at this site, which is found to be significant for certain wind directions (up to ±134 W m −2 ).
Q H from periods when the surface is dry is compared to Q H immediately after precipitation and measurements show that Q H is reduced by 27% when the surface is wet.This, along with observed relations between Q H and vegetation landcover fraction and spatial precipitation data, is used to estimate Q E .This new technique to model Q E results in an annual total within 1.4% of simultaneous independent eddy covariance measurements in the scintillometer source area.Modelled Q E and observed Q H and Q * are then used to determine monthly energy balance residual Q R .
The Q R is found to vary linearly with F CO2 and HDD.Annual Q F is calculated from the annual cumulative Q R and mean Q F for this location is 4.0 MJ m −2 day −1 , or 46.3 W m −2 .These figures are in line with other reported values of Q F in the literature.
Results from this work advance understanding of how spatial and temporal patterns of surface-atmosphere energy exchange are influenced by heterogeneous landcover, anthropogenic activity, and surface conditions (e.g.moisture availability) at 10 0 to 10 1 km length-scales.This information has potential to improve air pollution dispersion and transport models and enhance risk assessment for exposure to extreme temperatures in cities. Methodologically, spatial analysis techniques developed here represent a new path forward for eliciting additional information from spatial flux datasets in heterogeneous environments.Scintillometers are also shown to be a valuable tool for observing spatially integrated fluxes in urban areas at neighbourhood to city scales.
This will enable a broader range of urban morphologies to be explored through measurements.In particular, observations are needed of areas with tall buildings.Worldwide, such areas are increasing in spatial extent but are minimally understood from a local-scale energy balance perspective.From this, we may start to address how models should incorporate relevant physical processes related to surface-atmosphere energy exchange.

Figure 1 .Figure 2 .
Figure1.Scintillometer configuration in central London.Shaded contours are annual average combined source area weightings for all three scintillometers and for the eddy-covariance (EC) system.The scintillometer source areas are aggregated as the equal-weighted average of source area weights for each pixel.Grid reference is m (easting, northing) UTM zone 31N.Land cover dataset is fromLindberg and Grimmond (2011).

cFigure 3 .
Figure 3. (a) Annual quality control flag distribution, where 'Data available' are the data used for analysis, and (b) seasonal data coverage for day (K↓>1 W m −2 ) and night (K↓<1 W m −2 ).Horizontal black lines represent the range of data coverage.Section 2.3 gives a detailed description of quality control procedures and flag descriptions.
c 2016 The Authors.Quarterly Journal of the Royal Meteorological Society published by John Wiley & Sons Ltd on behalf of the Royal Meteorological Society.Q.J. R. Meteorol.Soc.(2017)

c
2016 The Authors.Quarterly Journal of the Royal Meteorological Society published by John Wiley & Sons Ltd on behalf of the Royal Meteorological Society.

Figure 4 .
Figure 4. Linear comparison statistics for (a) eddy covariance versus scintillometers, (b) Q H calculation methods (An, Db, FC), and (c) scintillometer paths (P1, P2, P3).The y-axis labels are listed in independent-dependent order for the linear comparison.Horizontal black lines represent range of data values.Details of each comparison are provided in section 3.1.Section 2.3 gives descriptions of Q H calculation methods and section 2.2, Figure 1 and Tables 1-2 give instrument configuration and site description.
(b)).c 2016 The Authors.Quarterly Journal of the Royal Meteorological Society published by John Wiley & Sons Ltd on behalf of the Royal Meteorological Society.

Figure 5 .Figure 6 .
Figure 5.Diurnal course of Q H for all three scintillometers during JJA.Boxplot shows data from all scintillometer paths combined and lines with symbols denote individual scintillometer path hourly median fluxes.The median hourly radiation and eddy-covariance (EC) Q H measured simultaneously at KCL are shown for reference.
c 2016 The Authors.Quarterly Journal of the Royal Meteorological Society published by John Wiley & Sons Ltd on behalf of the Royal Meteorological Society.

Figure 7 .
Figure 7.Diurnal course of seasonal mean Q H relative to the annual mean for all scintillometer pathlengths.DJF P1 data are not included due to extended periods of instrument maintenance (Figure3).

Figure 8 .
Figure8.Weekday:weekend Q H ratio of hourly medians for all seasons (box plots), winter (solid line) and summer (dashed line).P1 is excluded during winter due to insufficient data because of instrument maintenance (Figure3), and the dashed line is mean weekday:weekend ratio for JJA for all sensors.

Figure 9 .
Figure 9. (a) Frequency of observed Q H during JJA (all hours) at P2 in 20 • wind sector bins.Q H is binned in 100 W m −2 intervals, (b) Normalized Q H (Eq. (14)) in 20 • wind direction bins for all scintillometer sites, and (c) weekday:weekend Q H ratio of binned medians in 15 • bins.For (b)-(c), section 3.3 gives details on conditional sampling criteria to control for season, day of week, solar radiation input, and dry surface conditions.

Figure 10 .
Figure10.Normalized Q H (Eq. (14)) and source area landcover fraction of 20 • wind sector bins for (a) vegetation and (b) building for wet (white symbols) and dry (coloured symbols) surface conditions.Data points and linear statistics (a = slope, b = offset) are for wind sectors with n samples >10.Section 3.3 gives definitions of 'wet' and 'dry' surfaces and additional details on conditional sampling criteria to control for season, day of week, and solar radiation input.

Figure 12 .
Figure 12.Linear comparisons for monthly totals of (a) heating degree days (HDD) and Q R , (b) F CO2 measured by EC at KCL and Q R , and (c) HDD and F CO2 .Section 3.6 gives details on calculation of HDD and Q R .

Table 1 .
Scintillometer network configuration in central London.
c 2016 The Authors.Quarterly Journal of the Royal Meteorological Society published by John Wiley & Sons Ltd on behalf of the Royal Meteorological Society.

Table 2 .
Long-term source area-weighted surface landcover and morphological properties for scintillometer paths P1-P3 and the eddy-covariance (EC) method.
Quarterly Journal of the Royal Meteorological Society published by John Wiley & Sons Ltd on behalf of the Royal Meteorological Society.
c 2016 The Authors.Quarterly Journal of the Royal Meteorological Society published by John Wiley & Sons Ltd on behalf of the Royal Meteorological Society.