Efficiently modelling urban heat storage : an interface conduction scheme in an urban land surface model ( aTEB v 2 . 0 )

Intercomparison studies of models simulating the partitioning of energy over urban land surfaces have shown the heat storage term is often poorly represented. In this study, two implicit discrete schemes representing heat conduction through urban materials are compared. We show that a well-established method of representing conduction systematically underestimates the magnitude of heat storage compared with exact solutions of one-dimensional heat transfer. We propose an alternative method of similar complexity that is better able to match exact solutions at typically employed resolutions. The proposed interface conduction scheme is implemented in an urban land surface model and its impact assessed over a 15-month observation period for a site in Melbourne, Australia, resulting in improved overall model performance for a variety of common material parameter choices and aerodynamic heat transfer parameterisations. The proposed scheme has the potential to benefit land surface models where computational constraints require a high level of discretisation in time and space, for example at neighbourhood/city scales, and where realistic material properties are preferred, for example in studies investigating impacts of urban planning changes.


Introduction
The climate of cities differ from surrounding natural landscapes because urban structures change the terms of the surface energy balance (Oke, 1982).One important term is storage heat flux density ( Q S ): the net flow of heat per unit area into and out of materials.In highly urbanised ar-eas, Q S becomes the dominant term of the daytime urban energy balance (Cleugh and Grimmond, 2012).During the day, heat is absorbed by dense urban materials and at night is released back into the atmosphere.The altered behaviour of Q S fundamentally affects environmental processes such as atmospheric stability, the evolution of the boundary layer, convection, and pollution dispersion, and is key in establishing urban heat islands (Barlow, 2014).As climate affects human health, comfort and the energy use of buildings, many land surface models have been developed to simulate the impacts of different urban forms (Grimmond et al., 2009).However, the most recent urban model intercomparison project found that most participants significantly under-represented heat storage, with average midday Q S bias errors of −50 W m −2 (Best and Grimmond, 2014a).In a world of rapid urbanisation and changing global climate, researchers and planners are interested in efficient and accurate models that are able to investigate urban climate impacts.Improving heat storage representations at the neighbourhood-to city-scale will benefit those studies.This paper analyses the performance of a well-established method to represent heat storage in urban land surface models and compares it with a proposed alternative method.Although the proposed scheme is similar in complexity, the analysis undertaken here indicates it is more accurate at relevant resolutions and reduces overall energy flux errors of an urban land surface model compared with the established approach.
The paper is structured as follows.Section 1 introduces key issues in measuring and modelling urban heat storage.Section 2 describes the conduction schemes and evaluation methods.Section 3 presents scheme results for a simple ide-Published by Copernicus Publications on behalf of the European Geosciences Union.alised environment, and within an urban model and subjected to observational forcing.Section 4 discusses and concludes findings.

Measuring storage heat flux in cities
In simple environments it is possible to measure Q S directly through a network of heat flux plates (Nunez and Oke, 1977), but the variety of urban materials, orientations, skyview factors and internal building environments in a typical city make direct measurement at neighbourhood scales impractical.In order to practically account for the full diversity of materials in a heterogeneous urban landscape, the most commonly accepted method to measure Q S at the neighbourhood scale is by calculating heat storage as a residual of the urban energy balance (Roberts et al., 2006).Following Oke (1988) a full representation of the urban energy balance (in W m −2 ) is where inputs of net all-wave radiation flux (Q * ) and anthropogenic heat flux (Q F ) are balanced by sensible and latent turbulent heat fluxes (Q H and Q E ), net heat storage ( Q S ) and net heat advection flux ( Q A ). Using radiometers and eddy covariance techniques, Q * , Q H and Q E can be measured directly, and are considered representative of the neighbourhood if instruments are at a sufficient height above ground where effects of individual roughness elements are blended (Cleugh and Oke, 1986).Q A is often considered negligible and excluded if areas adjoining observation sites have similar urban characteristics (Roberts et al., 2006).Q F can be estimated through energy use and population density analysis (e.g.Sailor and Lu, 2004).Then, rearranging Eq. ( 1) with known terms on the right, residual heat storage is The residual approach has the inherent problem of accumulating all observational errors, which can be significant for turbulent fluxes (Wilson et al., 2002).Estimates of Q F and assumptions regarding Q A add to overall uncertainty.The daytime errors in Q S observations used in this study have previously been estimated at 23 %, doubling at night (Best and Grimmond, 2014b).Given the significant observational uncertainties, this study first evaluates conduction schemes with exact solutions to one-dimensional heat transfer, before comparing with observations of Q S calculated as a residual of the urban energy balance.

Simulating storage heat flux in cities
At scales of individual buildings, Q S of walls and roofs can be accurately modelled by solving the three-dimensional heat conduction equation at high resolution.However, simulating urban weather and climate over larger scales requires simplification in order to keep computation practical (Martilli, 2007).Some urban models replace heat transfer calculations with empirical functions of net radiation (Grimmond et al., 1991), or with idealised representations like the force-restore method (e.g.Porson et al., 2010).Another common method, used by most participants of the First International Urban Land Surface Model Comparison Project (PILPS-Urban), is to calculate a discretised version of the heat conduction equation (Grimmond et al., 2010).Best and Grimmond (2014b) found models using this method were generally better at simulating Q S , although all categories of models predicted Q S outside observational error in more than 50 % of daytime intervals.Overall, models partitioned too little energy into storage and too much into other fluxes, a result supported by a previous intercomparison at a different site (Grimmond et al., 2010;Best and Grimmond, 2015).Although other studies evaluating individual urban models have found good agreement between simulated and observed Q S (e.g.Roberts et al., 2006), PILPS-Urban showed that most models have difficulty predicting Q S at sites that had not previously been used to evaluate a model, where no parameter optimisation had been undertaken.This study suggests a method to improve Q S prediction using an urban model that is conceptually similar to many participants of PILPS-Urban.

Material thermal parameters
A potential weakness in evaluating urban model performance is the wide variety of parameters that could describe urban materials.Urban models are sensitive to material thermal parameter variation, particularly the magnitude and phase of storage heat flux and sensible heat flux (Oleson et al., 2008b).Selecting the best thermal parameters is not necessarily straightforward; as land surface models are a simplified representation of reality, model material parameters should only be viewed as abstract representations of observed physical quantities (Gupta et al., 1999).Researching and inputting realistic values based on local material parameters requires considerable effort, and can result in little improvement to performance, or even degrade it (Loridan and Grimmond, 2012).In PILPS-Urban, average model performance in Q S worsened as more site-specific (realistic) material thermal properties were provided to participants (Best and Grimmond, 2014b).As such, modellers sometimes use material parameters that do not match with the realities of the simulation site, but produce results that match well with observations (i.e.optimised parameters).This may pose problems for studies at new sites with different conditions, and studies that wish to ascertain impacts of changing urban materials.There is therefore interest in improving Q S prediction using more realistic material parameters.Additionally, modellers sometimes describe urban materials with a single set of thermal parameters (homogeneous materials), while others describe distinct layers with different thermal parameters (composite materials).In order to cover a range of modelling possibilities, we evaluate optimised, realistic, homogeneous and composite material parameters by drawing from five material dataset sources described in Appendix A.

Conduction representations
Three calculation methods of heat storage are compared: two discrete schemes and an exact solution.The two discrete schemes lump a material's heat capacitance at a temperature node and calculate solutions numerically at each time step.The exact method calculates continuous harmonic solutions to a periodic forcing.All three solve Fourier's law and the one-dimensional continuity equation where q is the conduction heat flux density (W m −2 ), λ the thermal conductivity (W m −1 K −1 ), T the temperature (K), d the depth (m), C is volumetric heat capacity (J m −3 K −1 ) and t is time (s).

Half-layer scheme
A common discretised approach is to locate the temperature node centrally within a homogeneous layer -a halflayer scheme (Fig. 1a).This approach is used in many urban land surface schemes, for example Town Energy Budget (TEB) (Masson, 2000), Single-Layer Urban Canopy Model (SLUCM) (Kusaka et al., 2001), Building Effect Parameterization (BEP) (Martilli et al., 2002), Community Land Model -Urban (CLMU) (Oleson et al., 2008a), Vegetated Urban Canopy Model (VUCM) (Lee and Park, 2008), and the Australian Town Energy Budget (aTEB) (Thatcher and Hurley, 2012).Models utilising this method vary in their spatial and temporal resolution, but typically resolve between 1 and 10 substrate temperature nodes, at between 5 and 60 min time steps (Grimmond et al., 2009(Grimmond et al., , 2010)).The half-layer method is based on well-established land surface models representation of thermal conduction through soil (e.g.Oleson et al., 2010), and is also used in multi-layer snow and sea ice models (e.g.West et al., 2016).
A discretised form of the conduction equation (Eq. 3) is where T m k is the temperature at the kth node at time step index m.Conduction between temperature nodes occurs through material layers with homogeneous and invariant thermal characteristics.Where two adjacent layers i and i + 1 have different depths d or conductivities λ, then the Layer 1 Layer 2 Layer 3 Figure 1.Conceptual diagram for two methods of discretising heat transfer through (three) homogeneous layers.(a) The half-layer scheme as commonly implemented in urban land surface models, and (b) the proposed interface scheme, which moves temperature nodes (T k ) from the centre of layers to the interfaces between them, with heat capacity (C i ) half of each adjacent layer.Discrete paths of conduction (q k,k+1 ) pass through layer resistances (R i ), while external/internal environment temperatures (T ext /T int ) and surface thermal resistances (R ext /R int ) control boundary fluxes (q ext /q int ).
half-layer scheme will represent an effective conductance (W m −2 K −1 ) between temperature nodes k and k + 1 as where resistance R i = d i /λ i (with i = k hereafter, i is dropped).Then, combining Eqs. ( 5) and ( 6) with the discretised form of the continuity Eq. ( 4) over time step length t, the general implicit formulation of the half layer scheme is Net heat storage flux density through to the next time step ( Q m+1

S
) is the sum of the change in energy stored in individual layers: Energy conservation is demonstrated by the equivalency of Eqs. ( 7) and (8) for layers 1 to n, as inner conduction fluxes (q m+1 1,2 . ..q m+1 n−1,n ) cancel, leaving heat storage flux density equal to outer conduction terms: www.geosci-model-dev.net/10/991/2017/Geosci.Model Dev., 10, 991-1007, 2017 where q m ext is the external admittance flux which drives transience in the system, q m int is the transmittance flux into the building, T m ext /T m int are external/internal environmental temperatures, and R ext /R int are external/internal surface thermal resistances.

Interface scheme
While the half-layer scheme lumps capacitance at the centre of layers, an alternative approach is to lump capacitance at the interface between layers (Fig. 1b).Since the paths of conduction between nodes are now completely within homogeneous layers, Eq. ( 6) simplifies to where i = k.The capacitance of the temperature node now takes half the value of the adjacent layers; in effect, we have swapped the need for an effective conductance for an effective capacitance.The general implicit formulation for the interface scheme is with outer temperature nodes represented by half the outer layer heat capacity only, i.e. at k = 1; C 0 d 0 = 0 and k = n + 1; C n+1 d n+1 = 0 for an n layer system.Total heat storage flux density for the interface scheme is Again, energy conservation is demonstrated through the equivalency of Eqs. ( 11) and ( 12), which after cancelling inner conduction terms leaves heat storage flux density equal to outer conduction terms:

Exact solution
Here we use the admittance procedure (Butcher, 2006;Davies, 1973), which calculates exact solutions to planar heat transfer through a series of homogeneous layers when subject to a steady sinusoidal forcing on one side, with a fixed temperature on the other.The international standard ISO 13786:2007 documents the method.The exact solution to the heat storage flux density Q S, exact for a composite wall is the net of outer conduction terms -the admittance flux density minus the transmittance flux density: where P is the period of sinusoidal forcing and H 12 and H 22 are components of a 2 × 2 complex-valued heat transfer matrix calculated via multiplication of individual heat transfer matrices over n homogeneous layers: where R ext and R int are external and internal surface thermal resistances and Z i ab a component of the ith layer heat transfer matrix, from outside in.Components of the ith heat transfer matrix are where j is the imaginary unit and is the periodic penetration depth of the ith layer.We use this method to calculate the exact time-varying Q S of a composite material over a forcing period.Periodic areal heat capacity (ISO, 2007a) is a useful measure of a composite materials ability to store heat over a sinusoidal cycle.It is a better measure than overall heat capacity or surface thermal admittance as it accounts for the periodic penetration depth of each material layer (for thick composites) as well as heat lost through transmittance (for thin composites).It can be calculated exactly as with units J m −2 K −1 .

Idealised evaluation methods
The half-layer and interface discrete schemes are compared with exact solutions of heat transfer.The external boundary is forced by T ext , a sinusoidal temperature variation of one degree over a 24 h period representing the combined effects of external environment temperature and incident radiation (or sol-air temperature): where T 0 = 290 K is the average external temperature.In this study, the internal boundary environment T int is fixed, equal to T 0 .Boundaries have fixed surface thermal resistances representing both convective and radiative transfer.Values are standard horizontal heat transfer rates set out in ISO6946, where R ext = 0.04 and R int = 0.13 m 2 K W −1 (ISO, 2007b).
For both discrete schemes, a linear system of equations describing the temperature evolution of each node is generated and solved by decomposition and back-substitution of a tridiagonal matrix (Thomas' algorithm).Discrete schemes are run for six forcing periods of 24 h, with the first five periods discarded as spin-up, and the last compared with the exact solution.For the exact solution, components of the heat transfer matrix are calculated numerically and harmonic solutions computed.Performance statistics are based those used in PILPS-Urban Phase 2, as described in Phase 1 (Grimmond et al., 2010).

Urban model evaluation methods
The discrete schemes are also compared within an urban land surface model (aTEB) forced by observations using the methodology of the PILPS-Urban Phase 2 (Grimmond et al., 2011).

Description of urban model: aTEB
The proposed interface scheme is implemented in the Australian Town Energy Budget (aTEB) urban land surface model (Thatcher and Hurley, 2012).aTEB was developed to act as the urban component of a regional or global climate model, so takes the highly efficient building-averaged approach where the generic urban unit is an infinite street canyon (Nunez and Oke, 1977).Canyon surfaces (walls, road, snow and vegetation) are connected to a bulk canyon air layer via an aerodynamic resistance network.Roofs and canyon air are then connected in parallel to the overlying atmosphere.
Although written from the ground up, aTEB is conceptually based on the influential Town Energy Budget (TEB) urban canopy model (Masson, 2000) with some modifications for Australian conditions.Modifications include the following: -In-canyon vegetation for suburban areas represented by a big-leaf model, adapted from Kowalczyk et al. (1994) but with a largely reduced set of prognostic variables.
-An air-conditioning component which pumps waste heat into canyons and prevents buildings acting as energy sinks during high-temperature periods, allowing energy closure at each time step.
-A two-wall canyon allowing radiative interactions between a sunlit and shaded walls, and a canyon airflow parameterisation with venting and recirculating regions, each integrated through 180 • for all possible street orientations, adapted from Harman et al. (2004a, b).
aTEB would be categorised as a "complex" urban model following the methodology of Grimmond et al. (2010Grimmond et al. ( , 2011)), primarily because of its canyon-based approach.Conceptually similar models include TEB (Masson, 2000), SLUCM (Kusaka et al., 2001) and CLMU (Oleson et al., 2008a).A significant differentiator amongst conceptually similar models is the parameterisation of heat exchange between canyon surfaces and turbulent air.By default, the SLUCM uses a form of the Jürges formula (1924), while TEB and CLMU use a form of Rowley et al. (1930).An alternative approach was developed by Harman et al. (2004b) where aerodynamic conductance is separately calculated for each canyon surface based on the airflow in different regions of the canyon for different canyon geometries.Urban canopy models that utilise forms of the Harman circulation scheme include the Single Column Reading Urban Model (SCRUM) (Harman and Belcher, 2006), the Met Office Reading Urban Surface Exchange Scheme (MORUSES) (Porson et al., 2010) and aTEB (Thatcher and Hurley, 2012).In order to assess how the proposed conduction scheme is affected by different aerodynamic heat transfer parameterisations, the Jürges, Rowley and Harman methods have been implemented in aTEB as described in Appendix C. Further details on aTEB is available in Thatcher and Hurley (2012) and Luhar et al. (2014).

Observational data
Observational data were obtained from flux tower measurements in a suburban site in Melbourne, Australia (Coutts et al., 2007a, b).Data include up-welling and down-welling long and shortwave radiation (K ↑ , K ↓ , L ↑ , L ↓ ), net all-wave radiation (Q * ), turbulent heat fluxes (Q H , Q E ), air temperature, pressure, wind, humidity and rainfall.Sampling rates were between 1 and 10 Hz, block averaged to 30 min intervals over 474.4 continuous days from 13 August 2003 to 28 November 2004.Measurements were taken at 40 m above ground, at a height where the effects of individual buildings were sufficiently blended so that measurements were considered representative of the neighbourhood.
The data used here are identical to that used in the First International Urban Land Surface Model Comparison Project www.geosci-model-dev.net/10/991/2017/Geosci.Model Dev., 10, 991-1007, 2017 (PILPS-Urban) Phase 2 (Grimmond et al., 2011), from which our evaluation methodology follows, that is: -Observed downwelling radiation, air temperature, pressure, wind, humidity and rainfall data were used as forcing data to run the urban model offline, i.e. without the need to be coupled to an atmospheric/earth system model.
-Observed upwelling radiation, turbulent heat and residual heat storage flux observations were compared with model output to evaluate the performance of the model.
-The initial 108.4 days of observation were treated as spin-up and excluded from analysis.
-The remaining 366 days were analysed, but if any flux were missing or gap filled in a time interval, all data in that interval were ignored, resulting in 8520 usable half-hour time intervals.
The site at Preston, Melbourne, is typical of low-tomedium density suburban housing in Australia, with detached one-to two-storey brick, timber and steel framed buildings, separated by roads, lawn and large trees.The site is classified in Best and Grimmond (2014b) as a local climate zone (LCZ) 6 (Stewart and Oke, 2012) or as an urban zone for energy exchange (UZE) medium density (Loridan andGrimmond, 2011, 2012).

Results
Conduction schemes are evaluated through two independent methods: (1) in a highly idealised environment to allow comparison to exact solutions to heat transfer and (2) within an urban land surface model and compared with observations.

High-resolution simulation
Discrete schemes are run at very high resolution (200 layers, each 1 mm deep) in order to test code validity.Table 1 shows the normalised standard deviation (SD) of high spatial resolution simulations over various time steps, where an SD of 1.0 means amplitudes of discrete and exact solutions match.At the higher time and space resolutions, both schemes converge towards the exact solution as expected.As the temporal resolution is decreased to 30 min time step amplitudes of Q S reduce in both schemes.

Realistic material parameters
We now test the performances of the discrete schemes at space and time resolutions more practical and typical of urban land surface models.The 20 wall and 12 roof configurations described in the CLMU database were tested at 30 min time steps over various levels of complexity.In Fig. 2, errors for 2-10 layer configurations are plotted individually (288 total), along with layer means.Of the two discrete schemes, the interface scheme has an average normalised SD closer to the exact solution for each layer up to ten (Fig. 2a).The mean normalised SD of Q S for the interface scheme begins above the exact solution at two layers, but decreases as resolution increases.For the half-layer scheme, SD begins below the exact solution and increases as resolution increases.For the interface scheme at 30 min time steps, mean spatial and temporal discretisation normalised SD errors almost cancel at four layers, but never cancel for the half-layer scheme.
Figure 2b shows the layer mean of the interface scheme mean absolute errors (MAE) are smaller than the half-layer scheme.The normalised MAE for the half-layer scheme monotonically decreases as the number of layers increases.Time and space discretisations both lead to negative bias in Q S , so MAE is minimised by increasing the number of layers indefinitely.On the other hand, the interface scheme average MAE declines until a minimum at four layers, where the positive bias of spatial discretisation approximately balances the negative bias of temporal discretisation.Using the same method for 1, 15, 30 and 60 min time steps, we find the optimal numbers of layers for the interface scheme are: 8, 6, 4 and 2, while the half-layer scheme only improves at higher resolutions (Appendix D).At higher resolutions, the interface scheme provides less benefit over the half-layer scheme, but is still able to simulate Q S more accurately for a range of realistic wall and roof assemblies (Fig. D1).
Having evaluated the performance of the two schemes over increasing material complexity, we evaluate performance for different composite thermal characteristics.Figure 3 shows the normalised SD and the MAE for each scheme at 30 min time steps for the 288 CLMU walls and roofs, plotted against the composite material's periodic areal heat capacity (Eq.20), herein called thermal mass.Materials with low thermal mass (e.g.lightweight framed walls or sheet metal) do not absorb and store as much heat over the diurnal cycle as materials with high thermal mass (e.g.concrete or brick).In contrast to discrete layers, thermal mass varies continuously, so the response of each scheme is represented as a locally weighted linear regression (LOWESS) (Cleveland, 1979).The thicker lines represent all 288 material configurations, while the thin dashed lines represent a subset of walls and roofs with four layers (the number of layers used later in the urban model analyses).3. Heat storage errors of the half-layer and the proposed interface schemes vs. exact solution over increasing periodic areal heat capacity (or thermal mass), with lightweight framed construction on the left, and heavyweight stone/ rubble construction on the right.For both (a) normalised standard deviation and (b) mean absolute error, the interface scheme locally weighted linear regression (LOWESS) outperforms the half-layer scheme overall (solid lines).Differences between schemes are more pronounced for walls/roofs with large thermal mass, and those represented with fewer layers (four-layers shown dashed).
Figure 3a shows the interface scheme LOWESS of normalised SD is closer to the exact solution for all values of thermal mass represented in the CLMU database.The halflayer scheme's LOWESS of normalised SD shows an increasing negative bias for larger thermal mass, while the interface scheme LOWESS is less steep.Figure 3b plots MAE (non-normalised) and shows schemes have greater difference in absolute errors for assemblies with higher thermal mass.Heat storage becomes a larger proportion of the energy balance in neighbourhoods with more heat capacity, so a scheme that is better able to represent Q S in those instances is beneficial.

Typically modelled parameters
We now evaluate the performance of the discrete schemes using wall and roof characteristics that are typical of previous urban modelling studies.Figure 4 displays the results of four walls, one realistic (SITE) and three optimised (WRF, UZE and aTEB).All walls are made up of four layers and later analysed in the aTEB urban model (Sect.3.2).The upper panels show Q S amplitude through a 24 h forcing period, while the lower panels show the error from the exact solution at each 30 min time step.In each case, the half-layer scheme under-simulates the diurnal amplitude of Q S , while the inwww.geosci-model-dev.net/10/991/2017/Geosci.Model Dev., 10, 991-1007, 2017  (Grimmond et al., 2011); (b) WRF: a homogeneous material which forms the default parameters for WRF land surface schemes (Chen et al., 2011); (c) UZE: a homogeneous material optimised to reduce WRF errors at 15 urban locations (Loridan and Grimmond, 2012); (d) an optimised composite used in aTEB (Thatcher and Hurley, 2012), based on the ECOCLIMAP database (Masson et al., 2003).For more information on material parameters, refer to Appendix A. Upper panels show storage heat flux density over a periodic cycle, lower panels show normalised error from exact solutions.The interface scheme reduces errors for these wall representations.
terface scheme matches more closely with the exact solution.
The mean absolute error (MAE) normalised by the absolute mean of the exact solution (Fig. 4 lower panel) is smaller for walls represented with the interface scheme.The same analysis is undertaken for roofs and shown in Fig. 5. Overall the interface scheme improves performance of these four-layer representations; however, the degree of improvement is smaller than the walls analysed above.Errors for both schemes are more pronounced in the UZE and aTEB roofs, which have greater total depths and therefore higher spatial discretisation errors (roof and wall characteristics listed in Appendix A).
As a consequence of the method of discretisation, the interface scheme has one additional temperature node over the half-layer scheme for any given number of layers (see Fig. 1).In Fig. 2, the interface average errors are lower than halflayer average errors for assemblies with n + 1 layers, implying the benefit of the interface scheme is broader than simply adding an additional node.A further test can be undertaken by dividing the first layer of a homogeneous material equally for the half-layer scheme only; then both schemes are represented by the same number of nodes, and the first node has the same heat capacity.Considering the homogeneous assemblies in Figs. 4 and 5 (WRF and UZE walls/roofs), in each case the five-layered (five-node) half-layer system has higher MAE than the corresponding four-layered (five-node) interface system.

Impact on heat storage
Figure 6 compares the net storage heat flux density ( Q S ) of observations and the urban model by calculating a mean flux for each hour of the day over the 12-month evaluation period.The method of determining observed Q S as the residual of the surface energy balance includes inherent uncertainty (see Sect. 1.1); however, the relatively long observation period used here reduces stochastic error.The upper panels show the mean hourly flux density over the diurnal cycle, the lower panels show error from observations.For each of the four datasets, the amplitude of the diurnal storage heat flux density increases when using the interface scheme compared to the half-layer scheme, similar to the behaviour seen in the exact analysis.As the half-layer method generally under-represents the magnitude of Q S , the interface scheme is better able to represent the observed storage heat flux over the full diurnal cycle.The interface scheme also improves simulated storage heat flux at the 25th and 75th quartiles, particularly during the day.
Half-hourly Q S performance statistics of aTEB using the Harman aerodynamic heat transfer method are listed in Table 2.For the four material databases, the interface scheme improves error metrics in most cases.These statistics are comparable to those reported in the PILPS-Urban Phase 2 intercomparison project, which used the same observational Upper panels show storage heat flux density over a periodic cycle, lower panels show normalised error from exact solutions.The interface scheme reduces errors for these roof representations, but improvements are less pronounced than for walls (Fig. 4).In each case, the interface scheme reduces errors by increasing heat storage amplitude.Only Harman aerodynamic heat transfer parameterisation shown (others qualitatively similar).When evaluation is limited to 3-month seasons, differences between schemes are qualitatively similar to annual results, although overall errors are higher in summer where radiation flux magnitudes are larger (not shown).The MAE for hourly averaged fluxes of a mean day shown here is different to MAE of each half-hour time step throughout the evaluation period (Table 2).
data (Grimmond et al., 2011).In that study, heat storage was the most poorly represented energy flux (Best and Grimmond, 2015), where the best Q S RMSE of any model in the final stage was 53 W m −2 , and the mean of participating models was 65 W m −2 .The interface scheme in general improves the urban simulation performance in Q S for these four datasets.However, an improvement in the performance in Q S may be offset by reduced performance in another flux of the energy balance (Loridan et al., 2010;Loridan and Grimmond, 2012), and so assessing the impact of the proposed conduction scheme on other energetic fluxes is necessary.

Impact on other fluxes
In Fig. 7, a Taylor diagram (Taylor, 2001) extends the statistical evaluation of the two conduction schemes to include sensible and latent heat, and upwelling longwave radiation fluxes.Shortwave radiation flux is omitted because it is unaffected by the conduction schemes, and downwelling radiation fluxes are prescribed.A Taylor diagram is useful in determining whether the difference error (measured in centred, or bias corrected, root mean square error: cRMSE) occurs from a change in variance (standard deviation: SD) or a change in pattern correlation (Pearson's correlation coefficient: r).cRMSE (as defined in Taylor, 2001) and SD are normalised by the standard deviation of corresponding flux observations to easily compare fluxes with different value ranges and variances.The two conduction schemes are tested with the four material datasets (Appendix A), along with three different aerodynamic heat transfer methods (Appendix C) for a total of 24 simulations.We also plot the anonymous results of participants in the PILPS-urban Phase 2 intercomparison (Grimmond et al., 2011).Simulations using the half-layer scheme are represented by unfilled coloured markers and the interface scheme with filled markers, with the two conduction schemes connected by a line.
For net storage heat flux ( Q S ), the interface scheme improves cRMSE for all simulations.The reduced cRMSE is primarily from improvement in normalised SD.Overall best performance is achieved with the aTEB default material parameters, followed by UZE, WRF and SITE.For each material dataset, the Harman heat transfer method performs better than the Rowley or Jürges methods.The length of the connecting line between conduction schemes compared to the tight grouping of the three heat transfer methods indicates that the conduction scheme generally has greater impact on performance of Q S than the different aerodynamic heat transfer methods.The spread of the material datasets indicates that the model is sensitive to material thermal parameter choices.
For upwelling longwave radiation flux (L ↑ ), the interface scheme improves cRMSE for all simulations.Improvements in both normalised SD and r contribute to the lower cRMSE.Rowley and Harman methods perform better than the Jürges method; however, the model is more sensitive to material dataset and the conduction scheme choice.
For sensible heat flux (Q H ), the interface scheme improves cRMSE for all simulations.The reduced cRMSE results primarily from improved normalised SD, as correlation is generally flat or slightly degraded.The Jürges and Harman methods perform better than the Rowley method.For Q H , the difference between conduction schemes is less pronounced than for Q S and L ↑ .
For latent heat flux (Q E ), the impact of the interface scheme on cRMSE is mixed, but on average degrades performance.No clear pattern emerges from either SD or r, and all simulations are tightly grouped.
The mean change in performance from half-layer to interface schemes for the 12 sets of experiments are presented in Table 3, with improvements in bold.A paired, two-sided t test for significance of the mean change in metric surpasses the 95 % confidence interval in each case.

Discussion and conclusion
We evaluated the performance of two implicit discrete schemes that represent heat conduction through urban materials.The half-layer scheme is well established and widely used in land surface models for urban structures, soil, snow and ice.It lumps heat capacitance at nodes centred within discrete layers.We proposed an alternative scheme of similar complexity that lumps heat capacitance at the interface between layers.We used two independent methods to evaluate the schemes: comparison with exact solutions to heat transfer in an idealised environment, and comparison with long time-series observations for an urban site with heat storage calculated as a residual of the urban energy balance.
In the idealised evaluation, a series of multi-layered assemblies of various complexities were subjected to a sinusoidal temperature forcing representing the diurnal cycle.The half-layer scheme was found to systematically underestimate Q S magnitude compared to exact solutions of one-dimensional heat transfer, while the interface scheme Anonymous results of participants in the PILPS-urban Phase 2 intercomparison are also plotted (Grimmond et al., 2011).Overall the interface scheme improves storage heat, sensible heat and longwave radiation fluxes.Latent heat flux response is less clear, but on average the interface scheme degrades performance.See Table 3 for summary statistics.
Table 3.Average performance statistics: change from half-layer to interface scheme for 12 pairs of experiments (improvement in bold).The significance value (p) of a paired, two-sided t test for the null hypothesis is given in parentheses; in all cases the 95 % confidence level is reached.

Flux
better matched exact solutions.For the half-layer scheme, time and space discretisation errors were both negatively signed, leading to under-representation of Q S .For the in-terface scheme, discretisation errors were oppositely signed, reducing average error.An optimal combination of space and time discretisation can be found to approximately cancel bias www.geosci-model-dev.net/10/991/2017/Geosci.Model Dev., 10, 991-1007, 2017 errors; for 30 min time steps the interface optimum was four material layers.Overall, the interface scheme provided greatest benefit for simpler representations (fewer layers), and for materials with larger periodic areal heat capacity (higher thermal mass).
In the urban model evaluation, we assessed the impact of implementing the interface scheme on performance of an urban land surface model (aTEB) over a 15-month observation period for a site in Melbourne, Australia.We evaluated four material parameter datasets and three common aerodynamic heat transfer parameterisations.In the urban model, the interface scheme tended to increase the diurnal magnitude of Q S .In all simulations, heat storage was under-represented compared with observations, so increasing Q S magnitude improved performance.Each flux is linked via the urban surface energy balance and so other fluxes were affected by a change in heat storage.As more energy was partitioned into storage, less was available for the turbulent fluxes and urban skin surface temperatures were lower.As both Q H and L ↑ were over-represented, the interface scheme improved performance in these fluxes too.For Q E , the impact of the interface scheme was mixed, but on average degraded performance.Material thermal characteristics were not chosen to achieve the lowest flux errors, but to assess the impact of the alternative conduction scheme for typical simulation setups.Four previously published material datasets were evaluated: SITE, WRF, UZE and aTEB.The interface scheme improved performance for all datasets; however, overall errors were still large.The SITE material parameters, which were supposed to match most closely with realities of the site, performed the worst.The optimised material datasets with unrealistically large storage capacities performed better.Three common canyon surface aerodynamic heat transfer methods were also evaluated: Rowley, Jürges and Harman.The Harman method had lower Q S errors than equivalent experiments using Rowley or Jürges, but for other fluxes there was no clear best method.Overall, the urban model was more sensitive to material thermal dataset and conduction method than to aerodynamic heat transfer method.
By physical reasoning, the interface scheme increases storage available to the transient external environment by representing heat capacity at skin surfaces, resulting in larger diurnal amplitudes of Q S compared with the half-layer scheme.This on average improves performance for idealised wall and roof assemblies based on real-world construction and material properties, and where urban land surface models are under-representing heat storage magnitude.However, the interface scheme will degrade the performance of urban models if the magnitude of Q S is already well represented, which may be the case in models where thermal properties of urban materials are optimised to increase thermal mass, which have higher temporal and spatial resolutions, or which have an altogether different representation of urban geometry.
Other than affecting flux predictions, the interface scheme can provide structural benefits to urban land surface models.The skin temperatures of urban surfaces are used in balancing energy budgets and determining radiant and turbulent fluxes.The interface scheme calculates skin temperatures prognostically, while models using the half-layer schemes diagnose skin temperatures as an additional calculation, or assume the first layer bulk temperature is representative of the skin temperature.For aTEB, moving from a half-layer to an interface conduction scheme avoided the additional calculations required to diagnose skin temperatures, and resulted in a 5 % reduction in average runtime for offline simulations.
In conclusion, the interface conduction scheme has the potential to benefit urban land surface models simulating environmental phenomena at scales that require a high level of discretisation in time and/or space for reasons of efficiency.Examples include numerical weather prediction, where many simulations are required in short timeframes, or climate studies that require simulation over long timescales.The interface scheme also improves performance in assemblies using realistic material thermal parameters, so may benefit large-scale studies investigating future impacts of urban design or climate mitigation measures.Results presented here are based on a single urban model with multiple configurations, and on a single observation site, so future work may extend evaluation to other sites and other urban models.aTEB (optimised, composite): From the ECOCLIMAP database (Masson et al., 2003) on which TEB (Masson, 2000) defaults are based, but with increased layer depths per Thatcher and Hurley (2012).The walls and roofs are not representative of typical building methods in Australia, but nonetheless give reasonable results for generic Australian cities.The roof is a layered composite with thermal conductivity and heat capacity of dense concrete, aerated concrete and insulation; the walls a composite of concrete and insulation and the road/soil of asphalt and dry soil.
Figure D1.Mean heat storage errors of the half-layer and the proposed interface schemes vs exact solutions over 2-10 layer representations and for 1, 15, 30 and 60 min time steps (per Fig. 2).For both (a) normalised standard deviation and (b) normalised mean absolute error, the interface scheme outperforms the half-layer scheme on average for the 288 representations of walls/roofs based on the CLMU database (Jackson et al., 2010).Optimum number of layers for each time step is identified with a star.The margin of improvement of the interface scheme is greatest at low time and space resolutions.
In Fig. D1 the analysis of Sect.3.1.2is repeated for 1, 15, 30 and 60 min time steps.For the half-layer scheme, the optimum combination of time step length and layer number is achieved by increasing the number of layers to 10+ for each time step length.For the interface scheme, an optimal number of layers can be achieved because space and time discretisation biases for storage heat flux are oppositely signed, resulting in normalised standard deviations closer to 1 (Fig. D1a) and lower mean absolute errors (Fig. D1b).In higher-resolution representations, the interface scheme advantage is relatively small.

Figure 4 .
Figure4.Heat storage response of half-layer and proposed interface schemes vs exact solutions for four wall types: (a) SITE: a composite of realistic materials at the Melbourne observation site(Grimmond et al., 2011); (b) WRF: a homogeneous material which forms the default parameters for WRF land surface schemes(Chen et al., 2011); (c) UZE: a homogeneous material optimised to reduce WRF errors at 15 urban locations(Loridan and Grimmond, 2012); (d) an optimised composite used in aTEB(Thatcher and Hurley, 2012), based on the ECOCLIMAP database(Masson et al., 2003).For more information on material parameters, refer to Appendix A. Upper panels show storage heat flux density over a periodic cycle, lower panels show normalised error from exact solutions.The interface scheme reduces errors for these wall representations.

Figure 5 .
Figure 5. Heat storage response of half-layer and proposed interface schemes vs exact solutions.Material datasets per Fig. 4 and Appendix A.Upper panels show storage heat flux density over a periodic cycle, lower panels show normalised error from exact solutions.The interface scheme reduces errors for these roof representations, but improvements are less pronounced than for walls (Fig.4).

∆QSFigure 6 .
Figure6.Hourly averaged diurnal heat storage response of an urban land surface model to four wall/roof datasets over 12 months.Material datasets per Fig.4and Appendix A. In each case, the interface scheme reduces errors by increasing heat storage amplitude.Only Harman aerodynamic heat transfer parameterisation shown (others qualitatively similar).When evaluation is limited to 3-month seasons, differences between schemes are qualitatively similar to annual results, although overall errors are higher in summer where radiation flux magnitudes are larger (not shown).The MAE for hourly averaged fluxes of a mean day shown here is different to MAE of each half-hour time step throughout the evaluation period (Table2).

Figure 7 .
Figure 7.A Taylor diagram of energy fluxes for each simulation with interface scheme (filled markers) and half-layer scheme (unfilled markers) connected by a line.Change toward a normalised standard deviation of 1 indicates improved variance (radial distance) and toward the bottom axis indicates improved pattern correlation (radial angle).Where both are improved, the centred root mean square (cRMSE) error decreases towards zero (interior arcs).Different colours indicate different material parameter datasets (SITE, WRF, UZE, aTEB; see Appendix A).Different marker shapes indicate different aerodynamic heat transfer parameterisations (Rowley, Harman, Jüges; see Appendix C).Anonymous results of participants in the PILPS-urban Phase 2 intercomparison are also plotted(Grimmond et al., 2011).Overall the interface scheme improves storage heat, sensible heat and longwave radiation fluxes.Latent heat flux response is less clear, but on average the interface scheme degrades performance.See Table3for summary statistics.

Table 1 .
Normalised standard deviation (SD) of Q S for highresolution (1 mm deep) simulations: aTEB wall test case.
(Jackson et al., 2010) the half-layer and the proposed interface schemes vs. exact solutions over various levels of complexity.Thirty-two roof and wall types from the CLMU database(Jackson et al., 2010)were modelled with nine levels of complexity, for a total of 288 representations (dots).Mean of layer errors indicate scheme performance at different levels of complexity (lines).For both (a) normalised standard deviation and (b) normalised mean absolute error, the interface scheme on average outperforms the half-layer scheme.Mean improvement diminishes as the number of represented layers increases.

Table 2 .
aTEB half-hourly performance statistics for Q S : different conduction schemes (improvement in bold).