WRF (v4.0)–SUEWS (v2018c) coupled system: development, evaluation and application

The process of coupling the Surface Urban Energy and Water Scheme (SUEWS) into the Weather Research and Forecasting (WRF) model is presented, including pre-processing of model parameters to represent spatial variability in surface characteristics. Fluxes and mixed-layer height observations in the southern UK are used to evaluate a 2-week period in each season. Mean absolute errors, based on all periods, are smaller in residential Swindon than central London for turbulent sensible and latent heat fluxes (Q_H, Q_E) with greater skill on clear-sky days on both sites (for incoming and outgoing short- and long-wave radiation, Q_H and Q_E). Clear-sky seasonality is seen in the model performance: there is better absolute skill for Q_H and Q_E in autumn and winter, when there is a higher frequency of clear-sky days, than in spring and summer. As the WRF-modelled incoming short-wave radiation has large errors, we apply a bulk transmissivity derived from local observations to reduce the incoming short-wave radiation input to the land surface scheme – this could correspond to increased presence of aerosols in cities. We use the coupled WRF–SUEWS system to investigate impacts of the anthropogenic heat flux emissions on boundary layer dynamics by comparing areas with contrasting human activities (central–commercial and residential areas) in Greater London – larger anthropogenic heat emissions not only elevate the mixed-layer heights but also lead to a warmer and drier near-surface atmosphere.

1 Introduction

Accurate prediction of urban–atmosphere interactions is one essential task of modern numerical weather prediction (NWP) models. There is increasing need to understand city-weather feedbacks and their impact on citizens and infrastructure to facilitate the delivery of integrated urban services (IUSs). IUSs span weather, climate, hydrometeorological and environmental processes and are related to many urban functions (e.g. urban planning, building design, transport/logistics operation, health, energy infrastructure and operations; Baklanov et al., 2018; Grimmond et al., 2020; Masson et al., 2020).

To improve such predictions, numerous efforts have been made to develop and enhance many urban land surface models (ULSMs; Grimmond et al., 2010a, b; Best and Grimmond, 2015), including the Single-Layer Urban Canopy Model (SLUCM; Kusaka et al., 2001), the Building Effect Parameterisation (BEP; Martilli et al., 2002), the Town Energy Balance model (TEB; Masson, 2000), and the Surface Urban Energy and Water Scheme (SUEWS; Järvi et al., 2011). Beyond resolving the transfer of energy, water and scalars at the land–atmosphere interface, the core tasks of ULSMs are to perform the following.

• Characterise the urban surface. The heterogeneous mix of materials and morphology varies from being dominated by built surfaces (i.e. buildings and paved areas) in city centres to having sparsely built fractions and more vegetation at more residential outskirts. Morphological variability is driven by the changing heights and spacings of buildings and trees across cities.

• Account for anthropogenic dynamics. As people's behaviour varies, it modifies emissions (e.g. energy, aerosols, water) on both regular (e.g. workweek and weekends) and irregular (e.g. major sports events such as Olympics, concerts, COVID19) patterns, modifying urban–atmosphere interactions. Thus, the city morphology or form remains relatively constant, but the functioning of the city varies with changing behaviour patterns.

• Capture the impact of urban–atmosphere interactions on the boundary layer. The urban boundary layer (UBL) is the lowest part of the atmosphere that is directly influenced by the presence of the city. The UBL is characterised by higher wind speeds and higher turbulent fluxes than the overlying free atmosphere. Additionally, it features elevated concentrations of pollutants and aerosols, as well as a warmer and more humid near-surface atmosphere.

SUEWS, a widely used and tested ULSM (Table 1), uses a mix of seven land cover types to characterise the surface materials. Anthropogenic heat, water and carbon emissions, with other features (e.g. snow clearing, irrigation), are used to capture behavioural dynamics impacts on urban–atmosphere interactions. Since its development, SUEWS has been regularly enhanced (e.g. Grimmond et al., 1986; Grimmond and Oke, 1991; Grimmond et al., 1991; Järvi et al., 2011, 2014, 2019; Offerle et al., 2003; Ward et al., 2016; Omidvar et al., 2022) and tested in a wide range of climates and cities worldwide (Table 1). Although operationally simple and scientifically robust, the full SUEWS model has primarily been used offline, preventing many urban–atmosphere feedbacks to be explored with the model. Coupling ULSMs (such as SUEWS) into larger-scale atmospheric models would better represent the land surfaces with more detailed physical processes and is thus expected to enhance the understanding of urban–atmosphere interactions (Vilà-Guerau de Arellano et al., 2023).

Omidvar et al. (2022)Tang et al. (2021)Lindberg et al. (2020)Kokkonen et al. (2019)Kokkonen et al. (2019)Sun and Grimmond (2019)Ao et al. (2018)Järvi et al. (2019)Kokkonen et al. (2018b)Kokkonen et al. (2018a)Lindberg et al. (2018)Ward et al. (2017)Demuzere et al. (2017)Järvi et al. (2017)Järvi et al. (2017)Rafael et al. (2017)Ward and Grimmond (2017)Alexander et al. (2016)Ao et al. (2016)Karsisto et al. (2015)Ward et al. (2016)Alexander et al. (2015)Onomura et al. (2015)Järvi et al. (2014)Järvi et al. (2011)Mitchell et al. (2008)

Table 1 Recent studies involving SUEWS have undertaken (1) development (D) of modules (M) and improvements of supporting tools (T) to coefficients (C), with (2) applications (A) where the model has been evaluated (E) or used to assess a scenario (S) outcome.

n/a: not applicable

Here, we couple SUEWS (v2018c; Sun et al., 2019) to the Weather Research and Forecasting (WRF) model (V4.0; Skamarock et al., 2019), an open-source frequently used NWP model. WRF provides the atmospheric forcing to SUEWS, and in turn WRF receives surface–atmosphere feedbacks for the city and the region. In this paper, we describe the structure and key physics of the coupled WRF–SUEWS system (Sect. 2), evaluate WRF–SUEWS at two UK sites (Sect. 3), and explore its application in modelling dynamics and impacts of anthropogenic heat emissions at the city scale (Sect. 4).

2 Development of the WRF–SUEWS coupled system

2.1 Physical interactions between WRF and SUEWS

The coupling between WRF and SUEWS occurs via the biophysical interactions between the land surface – with SUEWS introduced as a new land surface module option – and other physics modules in WRF (Fig. 1).

1. The radiation module provides radiative forcing variables, incoming short- (K_↓) and long-wave radiation (L_↓), to the land surface module. The land surface module returns outgoing short- and long-wave radiation (K_↑ and L_↑).

2. Atmospheric variables needed for SUEWS, including air temperature T_a, relative humidity RH, barometric pressure p_a and wind speed U, are supplied by the boundary layer (BL) module. These are influenced by turbulent transport (i.e. momentum τ, sensible heat Q_H and latent heat Q_E fluxes) from the land surface.

3. Precipitation P is generated by the microphysics and cumulus modules that parameterise precipitation-related processes at different scales.

Figure 1 Interactions of the five WRF physics (blue) schemes through processes (yellow) and the land surface module variables (purple). Notation defined in Sect. 2.2.

2.2 Key physics of SUEWS

SUEWS simulates both the energy balance (Oke, 2002),

$$\begin{array}{}\text{(1)}& Q^{*}+Q_F=Q_H+Q_E+\mathrm{\Delta }{Q}_S,\end{array}$$

and water balance (Grimmond et al., 1986),

$$\begin{array}{}\text{(2)}& P+I_{\mathrm{e}}=E+R+\mathrm{\Delta }S.\end{array}$$

The two are linked through the latent heat (Q_E) or evaporative E fluxes (by the latent heat of vaporisation). The water balance is driven by precipitation (P) and external water use (I_e). Whereas the surface energy balance is driven by the net all-wave radiation ($Q^{*}=K_{\downarrow }-K_{\uparrow }+L_{\downarrow }-L_{\uparrow }$, where K and L denote the short- and long-wave components, respectively, and arrows ↓ and ↑ in the subscript denote the incoming and outgoing directions, respectively) in all environments but additionally in cities, the human activities result in anthropogenic heat flux emissions (Q_F). The turbulent sensible (Q_H), latent (Q_E) and net storage heat flux (ΔQ_S); runoff (R); and change in water storage (ΔS) each have distinct responses that differ with land use and land cover. Traditionally, SUEWS has been mostly used for urban areas (Grimmond et al., 1986; Grimmond and Oke, 1991; Järvi et al., 2011; Ward et al., 2016), but for WRF–SUEWS it has been extended to non-urban contexts (Omidvar et al., 2022). Each model grid cell has up to seven land cover types (paved, buildings, deciduous trees, evergreen trees, grass/crops, bare soil and water) whose fractions and properties (e.g. height, albedo, leaf area index) can each vary between grid cells.

In WRF–SUEWS, Q_F is calculated using heating and cooling degree days (HDDs and CDDs) following the Sailor and Vasireddy (2006) approach:

$$\begin{array}{}\text{(3)}& Q_F={\mathit{\rho }}_{\mathrm{pop},t}\left[a_{F\mathrm{0}}+a_{F\mathrm{1}}\text{CDD}+a_{F\mathrm{2}}\text{HDD}\right],\end{array}$$

where a_F0,a_F1 and a_F2 are grid-specific coefficients. The grid population densities ρ_pop,t could be a daily mean (i.e. day and night) value (e.g. Ward et al., 2016) or capture the diurnal variations (e.g. Ward and Grimmond, 2017; Ao et al., 2018). Typically, for cities with strong commuting flows (e.g. London), Q_F and ρ_pop,t are larger in the central business districts (CBDs) during the day due to the commuting but are higher in residential areas at night. As using daily mean population density may bias Q_F and lose intra-daily variability (e.g. difference in large city centres between work-intensive and non-work periods), following other models (e.g. Allen et al., 2010) we divide the day into four periods (morning transition, day, afternoon transition and night). Day and night population densities are used in their respective periods, and their averages are used in both transition periods. The other parameter values can be derived from more detailed models that are not rapid enough for NWP (e.g. GQF, Iamarino et al., 2011; Gabey et al., 2018; DASH, Capel-Timms et al., 2020). ΔQ_S for each grid cell is calculated using the Objective Hysteresis Model (OHM) (Grimmond et al., 1991):

$$\begin{array}{}\text{(4)}& \mathrm{\Delta }{Q}_S=\sum _i{f}_i\left[a_{\mathrm{1},i}{Q}^{*}+a_{\mathrm{2},i}{\displaystyle \frac{\partial Q^{*}}{\partial t}}+a_{\mathrm{3},i}\right],\end{array}$$

where t is the time, f_i is the fraction of the area of each land cover type (i) and $a_{\mathrm{1}-\mathrm{3},i}$ values are material-related coefficients for each land cover type that can vary with the grid cell (cf. Sect. 3.1). This approach allows for much more rapid calculation of this flux than other methods. Q_E is introduced to the Penman–Monteith equation (Grimmond and Oke, 1991):

$$\begin{array}{}\text{(5)}& Q_E={\displaystyle \frac{s\left(Q^{*}+Q_F-\mathrm{\Delta }{Q}_S\right)+\mathit{\rho }{c}_{p}V/r_{\mathrm{a}}}{s+\mathit{\gamma }\left(\mathrm{1}+r_{\mathrm{s}}/r_{\mathrm{a}}\right)}},\end{array}$$

where s is the slope of saturation vapour pressure curve, ρ is the density of air, c_p is the specific heat capacity of air at constant pressure, V is the vapour pressure deficit, γ is the psychrometric “constant”, r_a is the aerodynamic resistance for heat or water vapour, and r_s is the surface or canopy resistance. Following Monteith (1965), assuming energy balance closure, Q_H is calculated as

$$\begin{array}{}\text{(6)}& Q_H=Q^{*}+Q_F-\mathrm{\Delta }{Q}_S-Q_E.\end{array}$$

With the latent heat of vaporisation (L_v) and Eq. (5) we obtain $E=\frac{Q_E}{L_{\mathrm{v}}}$ to link surface energy (Eq. 1) and water (Eq. 2) balance. R includes the runoff from individual surfaces, in channels and to groundwater. External water use (I_e) is estimated based on the automatic and/or manual irrigation or external application (e.g. street cleaning) as follows (Järvi et al., 2011):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{I}_{\mathrm{e}}& =A_{\mathrm{irr}}\left[f_{\mathrm{aut}}\left(b_{\mathrm{0},a}+b_{\mathrm{1},a}{T}_{\mathrm{d}}+b_{\mathrm{2},a}{t}_{\mathrm{r}}\right)\right.\\ & \left.+\left(\mathrm{1}-f_{\mathrm{aut}}\left(b_{\mathrm{0},m}+b_{\mathrm{1},m}{T}_{\mathrm{d}}+b_{\mathrm{2},m}{t}_{\mathrm{r}}\right)\right)\right],\end{array}\end{array}$$

where A_irr is the area irrigated, f_aut is the fraction of A_irr that is automatically irrigated, $b_{\mathrm{0}-\mathrm{2},a,m}$ values are site-specific coefficients, T_d is the daily mean temperature and t_r is the number of days since the rain. The net change in the water storage ΔS (e.g. in soil, in waterbodies, on the surface) is determined at each time step as the change in each surface water state compared to the previous time step.

The aerodynamic resistance (r_a) is calculated at first at the atmospheric level in WRF–SUEWS, where the wind speed (U) is determined (Fig. 1):

$$\begin{array}{}\text{(8)}& r_{\mathrm{a}}={\displaystyle \frac{\left[\ln \left(\frac{Z_{\mathrm{m}}-z_{\mathrm{d}}}{z_{\mathrm{0}m}}-{\mathit{\psi }}_{\mathrm{m}}\right)\right]\left[\ln \left(\frac{Z_{\mathrm{m}}-z_{\mathrm{d}}}{z_{\mathrm{0}\mathrm{v}}}-{\mathit{\psi }}_{\mathrm{v}}\right)\right]}{{\mathit{\kappa }}^{\mathrm{2}}u}},\end{array}$$

where z_d is the zero plane displacement height (m), z_0m (and z_0v) is the roughness length for the momentum (and heat/water vapour), u is wind speed at height Z_m, κ is the von Kármán constant (0.4) and ψ_m (and ψ_v) is the atmospheric stability functions for momentum (and water vapour). Stability is determined iteratively using the Obukhov length and initiated with a LUMPS-calculated (Grimmond and Oke, 2002) sensible heat flux taking the grid land cover fractions into account.

To compute the grid-integrated surface resistance (r_s), its inverse, the surface conductance (g_s), is used (Ward et al., 2016):

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{r}_{\mathrm{s}}^{-\mathrm{1}}& =g_{\mathrm{s}}=\sum _i\left(G_{\mathrm{PFT}}\left(g_{max,i}{f}_i{\displaystyle \frac{{\text{LAI}}_i}{{\text{LAI}}_{max,i}}}\right)\right.\\ & \left.\cdot g\left(K_{\downarrow }\right)g\left(\mathrm{\Delta }q\right)g\left(T_{\mathrm{a}}\right)g\left(\mathrm{\Delta }{\mathit{\theta }}_{\text{soil}}\right)\right).\end{array}\end{array}$$

The mix of vegetation within the grid is taken into account by considering each vegetation type i with land cover fraction f_i, the maximum conductance $g_{max,i}$, leaf area index (LAI_i), maximum LAI (LAI_max,i) and the surface conductance parameter (G_PFT) determined by plant functional type. Functions g(K_↓), g(Δq), g(T_a) and g(Δθ) are related to how the environmental variables – downwelling short-wave (K_↓) radiation, specific humidity deficit (Δq), air temperature (T_a) and soil moisture deficit (Δθ_soil ) – control the surface resistance. These functions have the following forms (Ward et al., 2016):

$$\begin{array}{}\text{(10)}& {\displaystyle }g\left(K_{\downarrow }\right)={\displaystyle \frac{\frac{K_{\downarrow }}{G_{\mathrm{K}}+K_{\downarrow }}}{\frac{K_{\downarrow ,\text{max}}}{G_{\mathrm{K}}+K_{\downarrow ,\text{max}}}}},\text{(11)}& {\displaystyle }g\left(\mathrm{\Delta }q\right)=G_{\mathrm{q},\text{ base }}+\left(\mathrm{1}-G_{\mathrm{q},\text{ base}}\right)G_{\mathrm{q},\text{ shape}}^{\mathrm{\Delta }q},\text{(12)}& {\displaystyle }g\left(T_{\mathrm{a}}\right)={\displaystyle \frac{\left(T_{\mathrm{a}}-T_{\mathrm{L}}\right){\left(T_{\mathrm{H}}-T_{\mathrm{a}}\right)}^{T_{\mathrm{c}}}}{\left(G_{\mathrm{T}}-T_{\mathrm{L}}\right){\left(T_{\mathrm{H}}-G_{\mathrm{T}}\right)}^{T_{\mathrm{c}}}}},\text{(13)}& {\displaystyle }g\left(\mathrm{\Delta }{\mathit{\theta }}_{\text{soil }}\right)={\displaystyle \frac{\mathrm{1}-\exp \left(G_{\mathit{\theta }}\left(\mathrm{\Delta }{\mathit{\theta }}_{\text{soil}}-\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{WP}}\right)\right)}{\mathrm{1}-\exp \left(-G_{\mathit{\theta }}\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{WP}}\right)}},\end{array}$$

where the G parameters are related to environmental controls indicated by subscripts K for solar radiation, q for the specific humidity deficit (“base” and “shape” for base value and curve shape, respectively), T_a for air temperature and θ for the soil moisture deficit. Table 2 gives the values used in the evaluation of the coupled WRF–SUEWS system (detailed in Sect. 3.2). $K_{\downarrow ,max}$ is the maximum incoming short-wave radiation (1200 W m⁻² used in this work); $T_{\mathrm{c}}=\frac{T_{\mathrm{H}}-G_T}{G_T-T_{\mathrm{L}}}$ with T_L and T_H being the lower and upper limits for switching off evaporation, respectively ($T_{\mathrm{L}}=-\mathrm{10}$ ^∘C and T_H=55 ^∘C); and Δθ_WP is the wilting point (120 mm).

LAI varies with growing degree days (GDDs) and senescence degree days (SDDs) via (Järvi et al., 2011, 2014)

$$\begin{array}{}\text{(14)}& {\text{LAI}}_{d,i}=\left\{\begin{array}{ll}& min\left({\text{LAI}}_{max,i},({\text{LAI}}_{d-\mathrm{1},i}{)}^{\mathrm{0.03}}\text{GDD}\right.\\ & \left.\times \mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}+{\text{LAI}}_{d-\mathrm{1},i}\right),\\ & T_{\text{Base,SDD}}<T_{\mathrm{a}}<T_{\text{Base,GDD}}\\ & max\left({\text{LAI}}_{min,i},({\text{LAI}}_{d-\mathrm{1},i}{)}^{\mathrm{0.03}}\text{SDD}\right.\\ & \left.\times \mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}+{\text{LAI}}_{d-\mathrm{1},i}\right),\\ & T_{\text{Base,GDD}}<T_{\mathrm{a}}<T_{\text{Base,SDD}}\end{array}\right.,\end{array}$$

where the previous-day (subscript d−1) LAI is used with the base temperature corresponding to the initiation of leaf-off (T_Base,SDD ) and leaf-on periods (T_Base,GDD). The model also requires ${\text{LAI}}_{max,i}$ and ${\text{LAI}}_{max,i}$ for each vegetation type i.

SUEWS accounts for the running water balance of the multiple surface types. The water amount on the canopy of each surface (C_i) (Grimmond et al., 1991) determines the surface resistance between dry and wet $\left(r_{\mathrm{s}}=\mathrm{0}\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}\right)$ by replacing r_s in Eq. (5) with r_ss (Shuttleworth, 1978):

$$\begin{array}{}\text{(15)}& r_{\mathrm{ss}}={\left[{\displaystyle \frac{W}{r_{\mathrm{b}}(s/\mathit{\gamma }+\mathrm{1})}}+{\displaystyle \frac{(\mathrm{1}-W)}{r_{\mathrm{s}}+r_{\mathrm{b}}(s/\mathit{\gamma }+\mathrm{1})}}\right]}^{-\mathrm{1}}-r_{\mathrm{b}}(s/\mathit{\gamma }+\mathrm{1}),\end{array}$$

where W is a function of the relative amount of water present on each surface to its water storage capacity S_i:

$$\begin{array}{}\text{(16)}& W=\left\{\begin{array}{ll}\mathrm{1}& C_i\ge S_i\\ \frac{K_{\mathrm{r}}-\mathrm{1}}{K_{\mathrm{r}}-S_i/C_i}& C_i<S_i.\end{array}\right.\end{array}$$

K_r depends on the aerodynamic and surface resistances,

$$\begin{array}{}\text{(17)}& K_{\mathrm{r}}={\displaystyle \frac{\left(r_{\mathrm{s}}/r_{\mathrm{a}}\right)/\left(r_{\mathrm{a}}-r_{\mathrm{b}}\right)}{r_{\mathrm{s}}+r_{\mathrm{b}}(s/\mathit{\gamma }+\mathrm{1})}},\end{array}$$

where r_b, the boundary layer resistance, is a function of friction velocity u_* (Shuttleworth, 1983):

$$\begin{array}{}\text{(18)}& r_{\mathrm{b}}=\mathrm{1.1}{u}_{*}^{-\mathrm{1}}+\mathrm{5.6}{u}_{*}^{\frac{\mathrm{1}}{\mathrm{3}}}.\end{array}$$

Equations (15)–(18) ensure that the surface resistance r_ss has a smooth transition from 0 s m⁻¹ (a completely wet surface) to r_s (a dry surface).

2.3 Major updates since SUEWS v2018c

This work presents a coupling framework and its evaluation using SUEWS v2018, ensuring consistency in internal physics with the offline version for the comparison later in Sect. 3.4. However we note that the coupling structure designed for WRF–SUEWS enables seamless upgrades to more recent SUEWS versions.

Current offline versions of SUEWS have options not in the coupled WRF–SUEWS system, including

• CO₂ fluxes for local-scale anthropogenic and biogenic urban–atmosphere exchanges (Järvi et al., 2019)

• roughness sub-layer profiles for the diagnosis of air temperature, humidity and wind speed within the roughness sub-layer (Theeuwes et al., 2019; Tang et al., 2021)

• 2-D radiation profiles for solar and thermal-infrared radiation for multi-layer urban canopies (Hogan, 2019)

• ESTM (elemental surface temperature method) for heat storage estimation using surface temperature and thermal properties (Lindberg et al., 2020).

2.4 Technical implementation of the WRF–SUEWS coupling

The following are considered in the design of coupled WRF–SUEWS system.

• Performance. Given coupling with file-based IO (input–output) exchanges has unacceptable computational performance, we use the SuMin module under the WRF framework (Figs. 2 and 3).

• Extendibility. As SUEWS is regularly enhanced (Table 1), it is desirable or even essential for the coupled system to use the full capacity of the standalone SUEWS.

• Sustainability. Given the vast community effort to build and improve sophisticated software systems, coupling should not be limited to one version. Instead a highly standardised coupling procedure is required to be sustainable (Meyer et al., 2020).

To address these, the coupled WRF–SUEWS system uses an adaptive intermediate layer SuMin (SUEWS in minimum mode) to link both models (Fig. 5). From SUEWS, SuMin calls the main SUEWS calculator to conduct all core SUEWS physics calculations. Whereas from WRF, SuMin is linked to the module_sf_suews via suews_1d as a complete land surface model that can be used by the WRF dynamics solver (i.e. ARW solver; Advanced Research WRF) via the surface driver (Fig. 5). By coupling SUEWS and WRF this way, fast prototyping of new functionalities is possible on the SUEWS side while maintaining a stable coupling to the more complex WRF. When new SUEWS features are available to be fully coupled, appropriate switches can be activated to incorporate them within the whole WRF system. This intermediate-layer-based approach allows for efficient communication between SUEWS and other models (e.g. SuPy, Sun and Grimmond, 2019) through an explicit, unified interface. Thus, SUEWS can be potentially coupled to other weather/climate modelling systems (e.g. OpenIFS, ECMWF, 2021).

Figure 2 WRF–SUEWS system consists of the new module_sf_suews added into WRF (blue) to interact with SUEWS via SuMin (green) using input files pre-processed by the SUEWS pre-processor (yellow). File-based input–output flow (dashed) and runtime calling logic (solid lines) are shown. suews_init is a subroutine to initialise the SUEWS module coupled into WRF, while wrfbdy.nc and wrfout.nc are standard WRF boundary condition and output files, respectively.

A Python-based WRF–SUEWS pre-processor system (WSPS; Fig. 5) formats the data to allow the additional parameters not in standard WRF input files (e.g. input variables in wrfinput.nc- and namelist-based configurations; cf. IO workflow in Fig. 3) to be incorporated, rather than modifying the WRF Pre-Processing System (WPS). The WSPS can be used for offline model spin-up runs to obtain the appropriate required initial conditions for WRF–SUEWS. The files prepared are

• i.

  wrfinput.nc, which is a modified version of WRF inputs with initial model states and other static properties, and

• ii.

  namelist.suews, which is the global configuration for the SUEWS model.

To generate these, four types of inputs are needed (Fig. 3).

• i1.

  Standard WRF input files for WPS. The geographic and meteorological data are processed by WPS to produce wrfinput.nc files for the model domains and provide the template for the WSPS.

• i2.

  Additional input files. The static SUEWS-specific properties (e.g. land cover, population density, building morphology) and optional files (e.g. suitable default parameters to be used when known values are unavailable) will precede the same information, if available, in the standard WRF input files within the coupled WRF–SUEWS system. Note that in the London context this is not required (see later in Sect. 3).

• i3.

  Standard SUEWS input files. These are files used by SuPy (Sun and Grimmond, 2019) for offline spin-up simulations to obtain appropriate model initial conditions (an example shown in Sect. 3.2). The SUEWS settings (e.g. physics options, population density profiles) are used to create the namelist.suews global settings.

• i4.

  Land cover reclassification settings. In namelist.suews the relations between land covers for WRF and SUEWS (Sect. 2.4) are prescribed.

The WSPS input files (Fig. 3) can have different spatial resolutions between files. The implemented netCDF processor obtains the static properties (i2) and initial condition (i3) and resamples them to the geospatial configuration (projection method, resolution and averaging strategy) of the base wrfinput.nc (i1) to produce the wrfinput.nc files for WRF–SUEWS. Subsequently, namelist.suews can be easily modified by hand without going through the WSPS if useful (e.g. to test different configurations for spin-up or change land cover mapping relations).

Figure 3 Workflow for the WRF–SUEWS pre-processing system (WSPS).

The seven land cover (LC) types can be assigned different parameter values (Table 2) per grid cell and/or can change with time. For example, the “grass” vegetation type can have varying parameters by season (e.g. rice–wheat rotation). The WSPS uses the namelist.suews global configuration file to translate the WRF land use (LU) data, e.g. IGBP-modified (International Geosphere-Biosphere Programme) MODIS (Moderate Resolution Imaging Spectroradiometer) with 20 LU classes or the USGS with 24 LU classes (NCEP, 2000) to the 7 LC classes (Fig. 4). Each SUEWS LC may combine fractions from multiple IGBP LU classes. Through this reclassification, WRF–SUEWS can use existing LU data for SUEWS simulations while allowing the other parameters to vary between grids. Note that a flexible number of WRF LUs (up to 100 in this release) can be specified to compose a SUEWS LC so that an extremely heterogeneous LU composition can be accounted for.

Figure 4 The WSPS can be used to reclassify the WRF–IGBP default MODIS 20-category land uses to SUEWS-specific land covers: specifically, the 20 MODIS land use categories (left) with 4 general categories (middle) are reclassified into 7 SUEWS land covers (right).

2.5 Bulk-transmissivity-based solar radiation correction

Incoming short-wave radiation (K_↓) is known to be overestimated by WRF because of unresolved clouds and/or aerosols (Jimenez et al., 2016; Lapo et al., 2017). However, if the forcing radiation is too large, the other surface fluxes and variables will be impacted. Thus, they should not compare well to observations, or alternatively if they do compare well, the variables are correct for the wrong reason. In urban areas, even on clear-sky days, there is often a large presence of aerosols that impact bulk transmissivities (e.g. Shanghai, Xu et al., 2011; Ao et al., 2018; Beijing, Dou et al., 2019; Sun et al., 2022; London, Ryder and Toumi, 2011; Kotthaus and Grimmond, 2018a; Warren et al., 2018). A bulk atmospheric transmissivity (τ) can be specified in the namelist.suews to partially correct the overestimation of K_↓ by WRF, which can be determined using K_↓ at both the top of atmosphere $\left(K_{\downarrow ,\mathrm{TOA}}\right)$ and the surface $\left(K_{\downarrow ,S}\right)$ are used (Oke, 2002):

$$\begin{array}{}\text{(19)}& \mathit{\tau }={\displaystyle \frac{K_{\downarrow ,s}}{K_{\downarrow ,\mathrm{TOA}}}}.\end{array}$$

As τ can vary seasonally (e.g. the cases in London and Swindon as shown in Table 4) we determine the median clear-sky difference (Δτ_c) between τ_WRF and τ_obs from the analysis of clear-sky days observations around the peak K_↓ (which occurs between 40 % and 60 % of the daylight hours). The K_↓ forcing (Fig. 1) for SUEWS in WRF ($K_{\downarrow ,\mathrm{W}-\mathrm{S}}$) is then corrected using the original one produced by WRF ($K_{\downarrow ,\mathrm{WRF}}$) as

$$\begin{array}{}\text{(20)}& K_{\downarrow ,\mathrm{W}-\mathrm{S}}=K_{\downarrow ,\mathrm{WRF}}-\left(\mathrm{1}-\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{c}}\right)K_{\downarrow ,\mathrm{TOA}}.\end{array}$$

Given the empirical nature of the parameter values, this correction can only be applied where observations are available. Here, we apply the correction to all time periods but separate the evaluation (Sect. 3.4.1) by sky conditions to assess effectiveness. Obviously, this simple correction is not a complete solution but rather an attempt to obtain more accurate K_↓ forcing for the coupled SUEWS and, hence, better surface feedbacks for the WRF atmospheric modules (Fig. 1).

3 Evaluation of WRF–SUEWS at two UK urban sites

3.1 Surface characteristics of evaluation sites

WRF–SUEWS is evaluated at the same two UK sites as in a previous SUEWS evaluation study (Ward et al., 2016; W16 hereafter) for consistency – these sites exhibit distinct urban characteristics:

• KCL. The King's College London Strand Campus (51^∘30^′ N, 0^∘07^′ W) is a dense central business district area in London (d03 of Fig. 5).

• SWD. Swindon (51^∘35^′ N, 1^∘48^′ W) is a residential area in the town of Swindon (d04 of Fig. 5).

WRF–SUEWS is set up with four nested model domains (Fig. 5) with grid spacing (domain and number of grids) being 9 km (d01, 100×100), 3 km (d02, 115×91) and 1 km (d03 and d04, 76×76).

Figure 5 Model domain configurations: (a) four simulation domains (d01–d04) and urban land cover (paved and buildings) fraction in d03 (1 km resolution) based on (b) the WRF default MODIS dataset and (c) updated information for Greater London from the URBANFLUXES project (Lindberg et al., 2020). The land cover information is accessible in Sun et al. (2023a).

As the default MODIS-based built (building and paved) fraction does not capture the surface heterogeneity within Greater London (Fig. 5b), we replace it with a high-resolution (i.e. 2.5 m) land cover map (Fig. 5c) derived from earth observation (EO) VHR (very high resolution) SPOT (Satellite pour l'Observation de la Terre) imagery (Mitraka et al., 2016; Marconcini et al., 2017) with more realistic surface information. This high-resolution dataset is processed using UMEP (Lindberg et al., 2018) to derive both land cover fractions for the seven SUEWS classes and other morphological parameters of roughness elements (e.g. building and vegetation heights, frontal area index; Table 2). The resulting dataset is upscaled to obtain 1 km resolution data for d03. As equivalent detailed land cover information is unavailable for d04, the land cover (plus building and vegetation height) for the single grid where Swindon site is located is modified based on values in Ward et al. (2016) (Table 2).

Table 2 Key parameters assigned in the four model domains (Fig. 5) vary with land cover that is impervious (paved: PAV, buildings: BDG) and pervious (deciduous trees and shrubs: DCT, evergreen trees: EVT; grass: GRA, crops: CRP; bare soil: BSO, water: WAT). Anthropogenic heat flux coefficients vary between weekday (WD) and weekend (WE) periods. Data sources are Chrysoulakis et al. (2018) (C18), Lindberg et al. (2020) (L19), Omidvar et al. (2022) (O22) and Ward et al. (2016) (W16). Inhabitant: inh.

^a $\mathrm{0.6}\le f_{\text{PAV}}+f_{\text{BDG}}\le \mathrm{1}$. ^b $\mathrm{0.16}\le f_{\text{PAV}}+f_{\text{BDG}}<\mathrm{0.6}$. ^c $f_{\text{PAV}}+f_{\text{BDG}}<\mathrm{0.16}$.

To help assign SUEWS parameters related to the surface characteristics, the land cover characteristics of the 1362 d03 grid cells within Greater London are analysed by plan area fractions of paved (f_PAV) and building (f_BDG) land covers (Fig. 6). The most common (N=171) LC grid combination (i.e. f_PAV−f_PAV 0.05 fraction bins) is predominately pervious (notably grass) with minimal impervious area (<0.05 for both paved f_PAV and buildings f_BDG). The second most frequent (N=112, f_PAV=0.15 and f_BDG=0.1) is also largely pervious. It is also noting that KCL and SWD (blue dots in Fig. 6) reside in densely built-up and moderately pervious domains, respectively, indicating the different nature in land cover composition. Because high-resolution property information is not readily available across the evaluation domains, the surface-related SUEWS parameters (e.g. albedo, emissivity, OHM coefficients) are simplified into three classes based on the paved and buildings fraction (f_PAV+BDG) from the gridded land cover (Table 2): (a) densely built areas (f_PAV+BDG>0.6) are assigned parameter values of KCL (W16), (b) suburban areas ($\mathrm{0.16}<f_{\text{PAV+BDG}}\le \mathrm{0.6}$) are assigned parameter values of SWD (W16) and (c) natural surfaces ($\mathrm{0}<f_{\text{PAV+BDG}}\le \mathrm{0.16}$) are assigned parameter values based on dominant vegetation (Omidvar et al., 2022). In doing so, we can utilise the available property data to accurately represent surface heterogeneity. Note that the WRF–SUEWS system allows for grid cell level surface characteristic parameters assignment (e.g. SUEWS simulation of Greater London by Lindberg et al., 2020).

Given the importance of population density to anthropogenic heat emissions (Eq. 3), output area day- and nighttime population data (UK ONS, 2013) are resampled to 1 km resolution for d03. For the grid point of SWD in d04, the Ward et al. (2016) values are used; otherwise, zero anthropogenic heat emission is set with population density being zero.

Figure 6 Frequency (colour, N) of land cover characteristics in d03 (Fig. 5) for Greater London (GL) (1362 grids of 1 km²). The individual grid cells are categorised first by the greatest land cover fraction with an impervious (IMP) split between paved surfaces (PAV) and buildings (BDG). Other fractions are deciduous trees (DCT), evergreen trees (EVT), grass (GRA), bare soil (BSV) and water (WAT). The blue dots indicate the cover around the KCL and SWD evaluation sites.

3.2 Model setup and spin-up

WRF–SUEWS is run with two-way nesting mode of 33 vertical levels (top at 5 kPa 11 layers in the boundary layer below 2000 m with lowest levels in d03 and d04 being ∼40 m a.g.l.; above ground level) for all four domains (Fig. 5). We note more vertical levels may be needed in detailed investigations of atmospheric features (e.g. temperature, precipitable water, etc.); here a moderate number of vertical levels are used as a balance between the computational cost and necessary representation of atmospheric profiles considering the focus of this work on the model development and evaluation of essential urban–atmosphere interactions. The atmospheric boundary conditions used are the 1^∘ × 1^∘ (latitude × longitude) National Centre for Environmental Prediction FNL (final) data (NCEP, 2000). The well-tested WRF “CONUS” (contiguous United States) physics suite (configuration since v3.9) is used with the land surface scheme changed to SUEWS (Table 3). The SUEWS physics schemes (Table 3, details provided in Sun et al., 2019) are selected for simplicity including using the building and tree heights with a rule-of-thumb method (Grimmond and Oke, 1999) for momentum roughness length and displacement height; additionally, snow and irrigation modules are turned off (following W16's KCL and SWD configuration). We use the WRF adaptive time step option to reduce the total run time while being numerically stable considering both the horizontal and vertical extent (Hutchinson, 2007). For us, the adopted time steps for domain 1 to 4 are around 72, 9, 3 and 3 s.

Thompson et al. (2008)Tiedtke (1989); Zhang et al. (2011)Iacono et al. (2008)Iacono et al. (2008)Janjić (1994)Janjić (1994)Grimmond et al. (1986); Grimmond and Oke (1991); Järvi et al. (2011); Ward et al. (2016)Kawai et al. (2009)Grimmond and Oke (1999)Dyer (1974); Van Ulden and Holtslag (1985); Högström (1988)Järvi et al. (2011)Loridan et al. (2011)Grimmond and Oke (1991)Järvi et al. (2014)

Table 3 Physics scheme in the WRF and SUEWS option tested for use in coupled simulations. The local (internal) option number is given in the column labelled “No.” RRTMG: Rapid Radiative Transfer Model for GCMs (general circulation models).

We evaluate WRF–SUEWS during the four seasons using 2-week periods in 2012 (Table 4). To generate appropriate initial conditions (i.e. model spin-up), we conduct offline SUEWS runs driven by observations collected at KCL and SWD (refer to the purple boxes in Fig. 1 and related notations in Sect. 2.2 for details about the atmospheric forcing variables as well as Table 2 for surface property settings) for 2012 until the soil moisture converges (<0.1 % difference in last time step between consecutive years): a period of 15 years is needed for KCL, while a period of 5 years is needed for SWD. For fully vegetated grids (i.e. f_DCT/EVT/GRA=1), the observations at SWD are used as forcing conditions to spin-up SUEWS for 4 years before convergence. The required initial states (i.e. soil moisture, leaf area index) are used for each WRF–SUEWS period (Table 4). For grid cells with f_PAV+BDG>0.6, the KCL-based initial states are prescribed to represent areas dominated by impervious surfaces, while SWD is used for the rest that are not completely vegetated. The appropriate complete pervious cover type values are assigned to the pervious cells based on their dominant land cover type.

Table 4 Time periods in 2012 used to evaluate WRF–SUEWS (Sect. 3.2) in London (KCL) and Swindon (SWD) with the observed air temperature (KCL at 49.6 m a.g.l., SWD at 10.6 m a.g.l.) and rainfall. Measurement details are given in Ward et al. (2013) and Kotthaus and Grimmond (2014a, b). Note that daylight saving time impacts all except for the January period; i.e. people's activities (e.g. work times) are 1 h earlier than UTC. Δτ_c is the median clear-sky transmissivity difference between τ_WRF and τ_obs. The clear-sky days are determined with a mean of τ_obs>0.8.

3.3 Evaluation data and metrics

The W16 evaluation of SUEWS (v2016a) uses 60 min radiation and turbulent fluxes observed at KCL and SWD (Ward et al., 2013; Kotthaus and Grimmond, 2014a, b). Both flux towers are located close to the centre (within a 200–300 m radius circle) of a 1 km² model grid cell. The source areas of the observed turbulent fluxes have a probable 50 % contribution from within ∼400 m of the flux tower at KCL (Kotthaus and Grimmond, 2014b) and a probable 80 % contribution from within ∼700 m at SWD (Ward et al., 2013). Here, we average the 30 min sample output from land surface scheme (i.e. SUEWS) to 60 min to compare with the observations.

The mixed-layer height (MLH), derived from continuous high-resolution (15 s and 10 m) attenuated backscatter observed with a Vaisala CL31 ceilometer at Marylebone Road (MR) in London (Kotthaus et al., 2016; Kotthaus and Grimmond, 2018a, b), is used to evaluate the model's ability to predict atmospheric boundary layer (ABL) dynamics. The MLH values have been compared to AMDAR, or Aircraft Meteorological Data Relay (the median difference between inversion heights and MLH is 346 m based on all time periods; for more evaluation results, refer to Kotthaus and Grimmond, 2018a). Various observations can be used to obtain the height of the boundary layer, but the results depend on the variable used (Kotthaus et al., 2023). For example, differences occur between using temperature inversion and MLH (e.g. at night, Kotthaus and Grimmond, 2018a) or between MLH and the turbulence-derived mixing height (MH; Kotthaus et al., 2018), the height where the vertical velocity variations falls below a threshold (Barlow et al., 2011; Halios and Barlow, 2017). On the other hand, when using the Mellor–Yamada–Janjić scheme, the WRF output PBLH (planetary boundary layer height) is derived from the height where turbulent kinetic energy falls below 0.2 m² s⁻² (Banks et al., 2016; Janjić, 1994).

Although the comparison of the aerosol-derived MLH from observations and the turbulence-based mixing height diagnosed from the model output (WRF PBL, hereafter referred to as WRF MH) may be affected to systematic differences (e.g. those associated with vertical resolution as suggested by Kotthaus et al., 2023), the comparable nature between MLH and MH enables the former to be a proxy to examine the latter modelled by WRF.

The evaluation metrics used with the number of data points (N) available from the model output (Y_mod) and observation (Y_obs) time series are the following.

1. Hit rate (HR).

   $$\begin{array}{}\text{(21)}& \mathrm{HR}={\displaystyle \frac{{\sum }_{j=\mathrm{1}}^{N}H\left({\mathit{\delta }}_{Y,j}-\left|Y_{\mathrm{mod},j}-Y_{\mathrm{obs},j}\right|\right)}{N}}\end{array}$$

   with Heaviside step function H defined by

   $$\begin{array}{}\text{(22)}& H\left(x\right)=\left\{\begin{array}{ll}\mathrm{0},& x<\mathrm{0}\\ \mathrm{1},& x\ge \mathrm{0}\end{array}\right.\end{array}$$

   and the threshold δ_Y,j being a value dependent on evaluation variable Y. We use the HR to evaluate the surface energy fluxes with δ_Yj=50 W m⁻² for radiative $\left(K_{\downarrow },K_{\uparrow },L_{\downarrow },L_{\uparrow },Q^{*}\right)$ and ${\mathit{\delta }}_{Y,j}=\mathrm{0.1}{Q}^{*}+\mathrm{50}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ for turbulent (Q_E and Q_H) fluxes (Hollinger and Richardson, 2005), respectively. If HR=0, it suggests none of model predictions are within the acceptable threshold set, while HR=1 indicates all fall into the acceptance range.

2. Mean absolute error (MAE).

   $$\begin{array}{}\text{(23)}& \mathrm{MAE}={\displaystyle \frac{{\sum }_{j=\mathrm{1}}^N\left|Y_{\text{mod }}-Y_{\mathrm{obs}}\right|}{N}}\end{array}$$

3. Mean bias error (MBE).

   $$\begin{array}{}\text{(24)}& \mathrm{MBE}={\displaystyle \frac{{\sum }_{j=\mathrm{1}}^N\left(Y_{\mathrm{mod}}-Y_{\mathrm{obs}}\right)}{N}}\end{array}$$

Both the MAE and MBE have units of the variable analysed (i.e. W⁻² for fluxes, m for MLH or MH) with an ideal value of 0 indicating perfect agreement with the observations. The MAE, unlike the root mean square error, treats all error equally (Willmott et al., 2017).

3.4 Evaluation results

3.4.1 Effect of bulk transmissivity correction on solar radiation

First, we evaluate WRF–SUEWS' skill at predicting incoming short-wave radiation (K_↓) as it is crucial to driving surface–atmosphere processes (Fig. 1). Given that fixing WRF's RRTMG radiation scheme (Iacono et al., 2008) tendency to overestimate K_↓ is beyond the scope of this study, we modify K_↓ for each grid to ensure that the land surface receives the appropriate energy to drive the land surface scheme (e.g. SUEWS) by accounting for the differences in bulk transmissivity on clear-sky days (Fig. 7; correction methodology in Sect. 2.5).

Generally, after the correction is applied the modified K_↓ agrees better with observations on clear-sky days (Fig. 8) with reduced overestimation. The improvement in the HR is minimal, but the MAE and MBE become smaller. This type of correction may cause underestimation of the peak K_↓ values in summer (Fig. 8c; cf. Fig. 8g), as the WRF overestimated transmissivity is larger in the afternoon (Fig. 7). Thus using a single bulk correction will have a diurnal bias in K_↓ from undercorrection (overcorrection) in the morning (afternoon). This is evident in the ascending trend of Δτ_Sim-Obs (cf. Fig. 7).

On cloudy days, the correction also improves the K_↓ performance. The MAE and MBE, as well as the HR, are enhanced, despite the correction parameter not being derived for cloudy conditions (Fig. 9). The MBE is reduced by more than 50 % for all simulation periods at both sites (except October 2012 at SWD) after correction. In general, we deem such correction necessary and effective for land surface modules in the WRF system (v4.0); hence we use K_↓-corrected simulation results throughout the following analyses.

Figure 7 Bulk transmissivity difference (model − observation) of 60 min median values (line) and the interquartile range (shading) during the daylight hours (normalised between sunrise (SR: 0) and sunset (SS: 1)) during four periods of the year (Table 4) at (a) KCL and (b) SWD. See Sect. 2.5 for correction details of Δτ.

Figure 8 Daytime 60 min incoming short-wave radiation (K_↓) fluxes during clear-sky days. Observed data are marked with median values (circles) and the interquartile range (IQR, vertical lines), while simulations by WRF–SUEWS are illustrated with median values (lines) and the interquartile range (shading). Both the original (solid) and corrected bulk transmissivity (dashed) are utilised across four 2-week periods (refer to Table 4) at the (a–d) KCL and (e–h) SWD sites. Scarce clear-sky periods occurred in April. For additional details on metrics and units, refer to Sect. 3.3.

Figure 9 Same as Fig. 8 but for cloudy days.

3.4.2 Online and offline simulated surface energy balance fluxes

Although our focus is on the online WRF–SUEWS system, it is useful to compare its performance to the offline (i.e. standalone SUEWS) version. As we force the latter directly with observations (refer to the “atmospheric forcing” variables in Fig. 1), removing potential forcing errors from the online system, we expect better performance and hence use this as a benchmark. The same simulation periods were used for the online and offline evaluation (Table 4). Given that the SUEWS v2018c kernel is almost the same as the v2016 kernel, the offline performance is very similar to the values reported by W16. The largest difference may be associated with an update to the OHM calculations, as v2016 uses the mean net all-wave radiation values of the 2 preceding hours, while v2018c uses the step-size-weighted average of two successive time steps.

Overall, the offline (upper row, Fig. 12) performance is better than online (lower row) at both sites for all fluxes. Although offline K_↓ should have no error, slight deterioration in performance occurs (Ward et al., 2016) because the high-resolution (e.g. 1, 5 min) forcing is not used but rather 60 min means are interpolated to 5 min and subsequently re-averaged to 60 min for evaluation. The offline K_↓ HR is greater than 0.93 for all four seasons at both sites, whereas the online K_↓ HR values are 0.25 to 0.75 (Figs. B1–B2) with HR<0.5 and MAE >60 W m⁻² in April and July (Fig. 12). The model performance in K_↑ is comparable at both sites (Fig. 12); the offline mode proves to be superior to the online mode – this is not surprising as it corresponds to the model's performance in K_↓.

The outgoing long-wave radiation L_↑ is modelled well (MAE <12.6 W m⁻²) both offline and online. But like K_↓, L_↓ is poorer for both cases (MAE online: 28.6; offline: 28.2 W m⁻²) and its diurnal patterns is quite captured by the model correctly – the offline mode lacks sufficient diurnal variability, while the online mode shows general underestimation (cf. Figs. 10, 11). Net all-wave radiation (Q^*) is impacted mostly by K_↓ during the day, making the WRF K_↓ correction important. The intra-annual range of the MAE for the online Q^* simulations is smaller for KCL (21–64 W m⁻²) than at SWD (37–72 W m⁻²), so the turbulent fluxes at SWD start with a potentially greater error.

Figure 10 Diurnal pattern of simulated (median: online – solid, offline – dashed; IQR: shading) and observed (median: circle, IQR: vertical line) incoming (a–d) short-wave (K_↓) and (e–h) long-wave (L_↓), outgoing (i–l) short- (K_↑) and (m–p) long-wave (L_↑), and (q–t) net all-wave (Q^*) radiation for four 2-week periods in different seasons (Table 4) at KCL. “N” indicates the number of hourly data points used in the analysis; other metrics are defined in Sect. 3.3.

Figure 11 As Fig. 10 but for the SWD site.

Clear-sky seasonality in the MBE for the turbulent heat fluxes (Q_E and Q_H) occurs at both sites (Fig. 12). At KCL there is a positive bias of varying magnitude (6 to 47 W m⁻²) for both fluxes, whereas at SWD Q_E is generally underestimated (−6 to −4 W m⁻²) and Q_H is overestimated (underestimated) by 25 W m⁻² (−31 W m⁻²) in January (July).

Overall, WRF–SUEWS better predicts Q_E than Q_H at both sites (${\mathrm{MAE}}_{Q_E}$ vs. ${\mathrm{MAE}}_{Q_H}$: 26 vs. 45 W m⁻² at KCL and 18 vs. 62 W m⁻² at SWD). The diurnal performance for the turbulent heat fluxes and Bowen ratio $\left(\mathit{\beta }=Q_H/Q_E\right)$ is similar for both offline and online runs for both KCL and SWD (Figs. 13, 14). The β indicates the turbulent heat fluxes are correctly partitioned, suggesting the model's robustness in turbulent heat flux partitioning even when there are variations in radiation accuracy (i.e. making the skill in simulating the absolute radiation fluxes less critical). The WRF–SUEWS daytime β agrees well with the observations at both KCL and SWD. However, when both fluxes are small (<10 W m⁻²) there are both larger observational errors (e.g. uncertainties due to nocturnal weak turbulence; Järvi et al., 2018; Mahrt et al., 2012) and ratios change rapidly. Under these conditions, the nocturnal β is overestimated (January at KCL; all seasons at SWD).

Figure 12 Model performance for (rows) radiative $\left(K_{\downarrow },K_{†},L_{\downarrow },L_{†}\right.$ and Q^*) and turbulent heat (Q_E and Q_H) fluxes for (rows) simulated online and offline modes at (columns) KCL and SWD assessed using (columns) three metrics (HR, MAE, MBE: colour – darker for poorer performance) for four periods (triangles). Triangles marked × indicate the model performance can be categorised unsatisfactory based on the one of following criteria: HR<0.5, MAE>40 W m⁻² and $\left|\mathrm{MBE}\right|>\mathrm{40}$ W m⁻². Metrics are defined in Sect. 3.3.

Figure 13 Diurnal pattern of simulated (median: online – solid, offline – dashed; IQR: shading) and observed (median: circle, IQR: vertical line) fluxes: (a–d) sensible (Q_H) and (e–h) latent (Q_E) heat fluxes and (i–l) Bowen ratio $\left(\mathit{\beta }=Q_H/Q_E\right)$ (hourly median) for (columns) four 2-week periods in different seasons (Table 4) at KCL. Metrics defined in Sect. 3.3.

Figure 14 Same as Fig. 13 but for SWD.

To set these results into context, it is useful to compare them to the results of other urban land surface models that have been evaluated in this region (e.g. SLUCM, Loridan et al., 2013; Tsiringakis et al., 2019; Best-1T and MORUSES with JULES, Hertwig et al., 2020). Loridan et al. (2013) (hereinafter L13) focussed on 3 June 2010 using KCL observations to evaluate the SLUCM (Kusaka et al., 2001) in WRF (Chen et al., 2011). Hertwig et al. (2020) (H20) evaluated two urban schemes with JULES (Best et al., 2011), Best-1T (Best et al., 2006) and MORUSES (Porson et al., 2010), over the 2011–2013 period at both the KCL and SWD sites. The evaluation strategies differ: L13 considers both online and offline performance but with different surface information; H20 is offline only but considers a range of configurations. For comparison we consider their more “advanced” configurations (L13: online, UZE – Urban Zones for Energy partitioning – for plan-area-index-based surface categorisation; H20: offline, a baseline configuration, CTRL-B, and a more sophisticated one, CTRL-M, with the more realistic anthropogenic forcing and detailed land cover information) using appropriate periods.

WRF–SLUCM performance at KCL for Q_E is better (WRF–SUEWS MAE=33 W m⁻²); note that the L13 study covers only 1 d (cf. our summer of 14 d). Similar to SUEWS, SLUCM's online mode performs slightly worse than its offline mode with respect to RMSE (online vs. offline), 82.8 vs. 82.6 W m⁻² for Q_H, while for Q_E, it is 16.4 vs. 14.2 W m⁻². At KCL, the annual offline MAE_SUEWS for Q_E is similar to CTRL-B $\left(\sim \mathrm{20}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\right)$ but larger than CTRL-M $\left(\sim \mathrm{15}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\right)$, whereas all three model configurations have similar performance at SWD $\left(\sim \mathrm{15}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\right)$. For Q_H, all three models are similar at KCL $\left(\mathrm{MAE}=\sim \mathrm{40}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\right)$, but at SWD, MAE_SUEWS is larger $\left(\sim \mathrm{30}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\right)$ than for the two H20 models $\left(\sim \mathrm{20}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\right)$.

Although not directly comparable, the SUEWS performance appears to be consistent with both offline and online performance of other urban land surface models in this area. Attribution of bias differences (e.g. different parameterisations, configurations, land cover information) is out of the scope of this study but is the focus of model comparison studies (e.g. Grimmond et al., 2010a, b).

3.4.3 Boundary layer depth

To assess the boundary layer depth the modelled MH (Sect. 3.3) is compared to the observed MLH at MR in London (Fig. 15). Generally, WRF–SUEWS underestimates daytime MH in all seasons except for winter (Fig. 15a) with MAE>300 m in the warmer periods (Fig. 15b, c). WRF–SUEWS slightly overestimates MH at night. The MBE varies between −195 and 252 m for the four periods. These values are smaller than WRF–SLUCM (with multiple PBL schemes, −288 to 539 m) simulations over Greater Paris evaluated using radiosonde observations (Kim et al., 2013). Similarly, evaluating WRF using eight different PBL schemes in Barcelona, Banks et al. (2015) found daytime MH to be underestimated (cf. elastic backscatter lidar), with the largest relative bias of −48 %. This is comparable to our summertime results (Fig. 15c).

Figure 15 As Fig. 13 but for MLH (observed) and MH (online, Sect. 3.3) at MR in London.

4 Application of WRF–SUEWS: impacts of anthropogenic heat (Q_F) on the atmospheric boundary layer

4.1 Modelled variability in anthropogenic heat emissions

The WRF–SUEWS model demonstrates satisfactory performance across various seasons and for two distinct urban areas (Sect. 3.1). Therefore, we employ it to examine the influence of Q_F on the atmospheric boundary layer – a unique characteristic of the coupled WRF–SUEWS system compared to the standalone SUEWS – in d03 (Fig. 5) during April 2012, before the Olympics disrupted typical patterns in July. As daylight saving time has begun by this time of the year, people's activities are shifted an hour earlier (e.g. typical workday is 08:00–16:00 UTC or 09:00–17:00 BST local time; British summer time). Peak Q_F emissions occur in central London (Fig. 16) during the daytime ($>\mathrm{300}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, Fig. 16c). By the late afternoon and through the night (20:00–05:00 UTC) the values in these areas (Fig. 16b) are smaller ($<\mathrm{130}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, Fig. 16d). The areas where large values occur differ between night and day.

At night the peaks occur across a larger area beyond central London, with some larger values towards the outskirts (Fig. 16b; cf. Fig. 16a). The nocturnal Q_F mean in d03 is smaller (17.1 W m⁻²) than the daytime (19.8 W m⁻²) but is more spatially consistent (Fig. 16d; cf. Fig. 16c).

These spatiotemporal patterns (Fig. 16a, b) are largely consistent with previous studies in London. Daily peak values differ between studies because of model grid cell size (Lindberg et al., 2013) and efforts over time to reduce carbon emissions and therefore energy use (Lindberg et al., 2013; Ward and Grimmond, 2017). Peak values (grid size) vary between 120 W m⁻² (3.2 km², Capel-Timms et al., 2020), $\sim \mathrm{150}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ (1 km², Hamilton et al., 2009; 1 km², Bohnenstengel et al., 2013) and 210 W m⁻² (3.2 km², Iamarino et al., 2011). Besides, the WRF–SUEWS-predicted average Q_F values for London (i.e. $\sim \mathrm{18}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$) are comparable with those estimated in other mega-cities (e.g. peak values of $\sim \mathrm{50}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ in the city of Osaka, Japan, by Narumi et al. (2009); annual average of $\sim \mathrm{20}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ in the Beijing–Tianjin–Hebei agglomeration by Feng et al. (2012); annual average of $\sim \mathrm{30}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ in Beijing by Yang et al. (2022)).

Figure 16 April (Table 4) weekday normalised anthropogenic heat flux in London (grey lines are boroughs) (a) during main working hours in the day (08:00–16:00 UTC) and (b) at night (20:00–05:00 UTC), with (c, d) the respective distributions of their actual values (in W m⁻²). Grids analysed in more detail are indicated (blue): commercial–business (C) and residential (R).

4.2 Feedbacks from anthropogenic heat emissions

To assess the feedback, we focus on two example grids with contrasting population density profiles:

• a central business district (CBD) area, within the City of Westminster (grid C, Fig. 16), with tall buildings (roughness length z₀=2.0 m, zero plane displacement z_d=13.9 m; calculated using the rule-of-thumb approach as in Grimmond and Oke, 1999), a low fraction of vegetation (f_VEG=13.1 %) and large daytime population density (763 inh ha⁻¹) but small nocturnal density (111 inh ha⁻¹)

• an inner-city residential area, within the Borough of Islington (grid R, Fig. 16), with z₀=0.8 m, z_d=5.6 m, f_VEG=43.3 % and nocturnal population density (170 inh ha⁻¹) slightly greater than the daytime density (150 inh ha⁻¹).

Workday time series reveal differences in Q_F timing and magnitude between the two grids (Fig. 17a). In grid C, Q_F is significantly higher in general, staying at a rather consistent level between the morning and evening peaks, whereas the flux is lower in grid R and shows a reduction in emission between the two peaks. However, the second peak in R is much later (R at ∼21:00 UTC, C at ∼17:00 UTC, Fig. 17a, red). Both grids have similar Q_E (Fig. 17a) despite the 30 % difference in vegetation fraction, attributable to the low LAI in April (Kotthaus and Grimmond, 2014a, b; Ward et al., 2013). The net radiation (Q^*, Fig. 17a) values are also similar. The larger Q_F in grid C enhances Q_H by ∼220 W m⁻² over grid R during the middle of the day on average, which contributes to an increase in midday MH of 100–200 m compared to grid R (700–900 m a.g.l.; Fig. 17b, c).

The contrasting Q_F-induced heating alters several atmospheric variables within the urban boundary layer: daytime near-surface air $(<\mathrm{100}\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.)$ in grid C is warmer (positive difference, ∼0.3 K, Fig. 17b) and drier (negative difference, $\sim \mathrm{0.2}\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}$, Fig. 17c) than in grid R. At night, the air temperature in grid C stays warmer at lower altitudes but with cooler and wetter air aloft with PBL (Fig. 17). These results are consistent with other WRF–SLUCM-based studies stating that urban heating leads to warmer and drier near-surface atmosphere (e.g. Zhang and Chen, 2014; Zhang et al., 2015). The Q_F impacts are of similar magnitude to those linked to urban greening (but with an inverse effect): for instance, with all roofs vegetated in Beijing, the near-surface atmosphere is cooled by ∼1 K but moistened by $\sim \mathrm{0.8}\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}$ during a heat wave period even when anthropogenic heat is accounted for (Sun et al., 2016); this suggests that urban greening may help mitigate some effects of anthropogenic heat emissions.

Figure 17 April (Table 4) weekday diurnal pattern in two grids (C and R, Fig. 16) of (a) mean surface energy fluxes and (b, c) median MH and median difference (colour, C − R) with the height of (b) potential temperature θ and (c) specific humidity q.

5 Concluding remarks

Through coupling the SUEWS urban land surface model to WRF, urban–atmosphere interactions can be explored more fully than when using the standalone SUEWS. The new Fortran subroutine SuMin interfaces SUEWS with the land surface driver of WRF. The WSPS pre-processor incorporates SUEWS-specific parameters into the wrfinput.nc and namelist.suews files for WRF–SUEWS simulations. The coupling is designed to permit sustainable development of SUEWS so that regular enhancements of SUEWS can be seamlessly incorporated into the coupled system.

Evaluation of the coupled WRF–SUEWS system is performed at two UK sites: dense central London (KCL) and a suburban–residential site in Swindon (SWD) across four seasons. It generally shows a good capacity of the coupled system to simulate the surface energy balance fluxes and mixing height in all periods. The performance compares well to other WRF studies in urban settings (Banks et al., 2015; Kim et al., 2013; Loridan et al., 2013). Better performance is found (i) at the suburban site SWD compared to the densely built-up KCL for the turbulent heat fluxes, (ii) during clear-sky conditions for radiative compared to turbulent heat fluxes and (iii) using a bulk atmospheric transmissivity for incoming short-wave radiation (compared to without).

WRF–SUEWS' capability allows for analyses of spatial and temporal variations over heterogeneous urban areas compared to existing WRF urban schemes (e.g. SLUCM, BEP, etc.) that can only resolve a limited number of urban classes. Critically, the SUEWS capability allows for dynamic feedbacks from human activities (e.g. heat, water, phenology, snow related). The influence of anthropogenic heat Q_F on the boundary layer in April prior to leaf growth in Greater London influences the sensible (more than the latent) heat fluxes. The larger Q_F and Q_H in central London are associated with warmer and drier air and deeper mixing heights during the day but not at night. The WRF–SUEWS evaluation should be expanded to other urban settings, time frames and synoptic conditions, with further applications being explored (e.g. Q_F impacts on urban–atmosphere interactions). Results suggest that the system has a great potential to help advance our understanding of the role of urban surface heterogeneity. As SUEWS is already integrated with many other models (e.g. building energy parameterisations and thermal comfort simulations, Table 1), WRF–SUEWS can help decision-makers involved in a wide range of integrated urban services to identify the spatiotemporal distribution of near-surface meteorology (e.g. heat-induced climate risks, Zong et al., 2022; urban ventilation issues, Zheng et al., 2022) across a city.