Evaluation of the SPARTACUS-Urban Radiation Model for Vertically Resolved Shortwave Radiation in Urban Areas

The heterogenous structure of urban environments impacts interactions with radiation, and the intensity of urban–atmosphere exchanges. Numerical weather prediction (NWP) often characterizes the urban structure with an infinite street canyon, which does not capture the three-dimensional urban morphology realistically. Here, the SPARTACUS (Speedy Algorithm for Radiative Transfer through Cloud Sides) approach to urban radiation (SPARTACUS-Urban), a multi-layer radiative transfer model designed to capture three-dimensional urban geometry for NWP, is evaluated with respect to the explicit Discrete Anisotropic Radiative Transfer (DART) model. Vertical profiles of shortwave fluxes and absorptions are evaluated across domains spanning regular arrays of cubes, to real cities (London and Indianapolis). The SPARTACUS-Urban model agrees well with the DART model (normalized bias and mean absolute errors < 5.5%) when its building distribution assumptions are fulfilled (i.e., buildings randomly distributed in the horizontal). For realistic geometry, including real-world building distributions and pitched roofs, SPARTACUS-Urban underestimates the effective albedo (< 6%) and ground absorption (< 16%), and overestimates wall-plus-roof absorption (< 15%), with errors increasing with solar zenith angle. Replacing the single-exponential fit of the distribution of building separations with a two-exponential function improves flux predictions for real-world geometry by up to half. Overall, SPARTACUS-Urban predicts shortwave fluxes accurately for a range of geometries (cf. DART). Comparison with the commonly used single-layer infinite street canyon approach finds SPARTACUS-Urban has an improved performance for randomly distributed and real-world geometries. This suggests using SPARTACUS-Urban would benefit weather and climate models with multi-layer urban energy balance models, as it allows more realistic urban form and vertically resolved absorption rates, without large increases in computational cost or data inputs. Supplementary Information The online version contains supplementary material available at 10.1007/s10546-022-00706-9.


Introduction
Given their high concentrations of both people and infrastructure, cities are places of high vulnerability to variations in weather, climate, and air quality (Baklanov et al. 2018). Currently, limited-area numerical-weather-prediction (NWP) models have spatial resolutions such that cities span multiple model grid boxes (e.g., Hagelin et al. (2017)). In addition, global NWP models generally have poor representation of urban structure and energy exchanges, for example the European Centre for Medium-Range Weather Forecasts IFS (Integrated Forecasting System) model with 9-18 km resolution, with no urban model (Hogan et al. 2017). As NWP spatial resolution increases, smaller-scale processes will need to be resolved. Alongside this, global population and urban land cover are expected to increase, bringing about a need for greater understanding of energy exchanges between the surface and the atmosphere (Loridan and Grimmond 2012).
Cities have complex three-dimensional structures, with varying building heights, densities, materials, arrangements, shapes, and surroundings. These affect the radiative exchanges and heat storage within the urban surface (Grimmond et al. 2010;Yang and Li 2015;Ao et al. 2016), e.g., altering the effective shortwave albedo due to multiple reflections in the street canyon (Aida and Gotoh 1982). Given that shortwave radiation is the most important contribution to the surface energy balance, it is vital to understand how it is absorbed within urban areas (Fortuniak 2008).
Models have been developed to account for the three-dimensional nature of urban surfaces. A common approach is to simplify the urban form as a canyon of infinite length between buildings that are of equal height with a fixed height-to-width (H/W ) ratio (Nunez and Oke 1977). This urban-canyon approach is fast enough for NWP, and has been applied to urban radiation specifically (e.g., Aida 1982;Arnfield 1982aArnfield , 1988Harman et al. 2004) and other energy balance fluxes (e.g., Masson 2000;Kusaka et al. 2001a;Lee and Park 2008). This approach subdivides the canyon into three facets: walls, roof, and ground. Advancements of this approach include models subdividing facets further into sunlit and shaded, with varying canyon orientation (e.g., Oleson et al. 2008a, b), accounting for different building heights (Martilli et al. 2002), and the interactions between neighbouring canyons (Schubert et al. 2012). Some models have added vegetation, both at ground level and in the vertical plane (e.g., street trees) (Lemonsu et al. 2012;Krayenhoff et al. 2014;Redon et al. 2017).
These improvements in the vertical structure of urban form have led to improvements in the prediction of shortwave fluxes onto roofs (Schubert et al. 2012). Despite this, many models still make the unrealistic assumption of an infinitely long urban canyon. This leads to models neglecting key features of the urban form, such as intersections, building height variations, courtyards, and clusters of buildings. Ignoring these features impacts building shadowing, radiation trapping between buildings, and increased penetration of shortwave radiation to the surface in open areas such as parking areas, hence impacting the overall energy balance.
Building-resolving models, with details of each individual building and facet, are suitable for microscale research applications but not NWP, given their high data and computational demands (e.g., Krayenhoff and Voogt 2007;Krayenhoff et al. 2015;Resler et al. 2017). These models simulate radiative interactions between individual buildings, requiring detailed three-dimensional (3D) geometry and material data, which are hard to obtain for large areas (Ghandehari et al. 2018;Masson et al. 2020). Gastellu-Etchegorry et al. (2012) suggested a key application for complex building-resolving models is both to calibrate and evaluate simpler radiative transfer models (e.g., suitable for NWP). However, very few urban radiative transfer models have been evaluated against these models. One exception is the evaluation of shortwave radiation in the Town Energy Balance (TEB) model against the SOLENE explicit radiative transfer model, considering vertical vegetation (Redon et al. 2017), although an infinite street canyon is used in both the TEB and SOLENE models. Other urban radiative transfer models suitable for NWP could (or should) be calibrated and/or evaluated using explicit 3D models.
Here, we evaluate the shortwave performance of the SPARTACUS (Speedy Algorithm for Radiative Transfer through Cloud Sides) radiative transfer model for urban areas, SPARTACUS-Urban (Hogan 2019a), with respect to the more detailed explicit 3D Discrete Anisotropic Radiative Transfer (DART) model (Gastellu-Etchegorry 2008). The SPARTACUS-Urban model resolves the vertical structure of the urban canopy and exploits the recent finding that wall-to-wall separation distances in an urban area fit an exponential distribution well (Hogan 2019b). It can account for atmospheric absorption, emission, and scattering in the urban environment, rather than assuming a vacuum, as is done by most urban radiation models (e.g., Masson 2000;Harman et al. 2004), while being fast enough for use in NWP. Although it can represent vertical profiles of vegetation, this capability is not evaluated here. We examine SPARTACUS-Urban's ability to predict the vertical profile of the clear-air downwelling and upwelling fluxes, and the absorption into walls and roofs, across a range of urban forms: from simple cuboid 'buildings' to highly realistic structures.
The methods include a description of the SPARTACUS-Urban (Sect. 2) and DART (Sect. 3) models, and of the evaluation techniques (Sect. 4). After an investigation of the underpinning assumptions that SPARTACUS-Urban makes about the urban form (Sect. 5), the results of the evaluation are presented (Sect. 6). Finally, a comparison of SPARTACUS-Urban with DART is made with respect to the Harman et al. (2004) method for radiation within an urban street canyon (Sect. 7).

Description of the SPARTACUS-Urban Model
The SPARTACUS-Urban model (Hogan 2019a) uses an approach that originates from the SPARTACUS model for simulating 3D radiative transfer in complex cloud fields (Hogan and Shonk 2013;Hogan et al. 2016) and SPARTACUS-Vegetation for 3D interaction of radiation in forest-type vegetation (Hogan et al. 2018). These algorithms share a common mathematical approach for treating radiative transfer in the presence of objects that are randomly distributed in the horizontal. The SPARTACUS-Surface open-source software package (Hogan 2021) combines the capabilities of SPARTACUS-Urban and SPARTACUS-Vegetation, but since vegetation is not considered here we refer to the algorithm as SPARTACUS-Urban.
The SPARTACUS-Urban model is underpinned by the one-dimensional discrete ordinate method (Stamnes et al. 1988), which assumes that diffuse radiation travels in 2N streams of different elevations, with N streams per hemisphere. As N increases, the radiation field is described more accurately, but with increased computational cost. Here 16 streams (N = 8) are used. The SPARTACUS-Urban model splits a scene (defined here as any combination of building geometry, solar zenith angle, and albedo) vertically by height, z, into n horizontal layers above an assumed flat ground level. Each layer is split horizontally into 'regions' of clear sky, vegetation, or buildings. As with other urban models, SPARTACUS-Urban computes the interaction of radiation with three urban facets (wall, roof, ground) and optionally vegetation.
Through the process of this work, we have found that two modifications to SPARTACUS-Urban are needed. Section 2.1 describes how to treat the more common occurrence of large open spaces such as parking areas and parks than predicted by the exponential model of urban geometry used by SPARTACUS-Urban. Section 2.2 describes a correction to account for fine structure in the building perimeters.

Modification to Treat Non-Exponential Building Separations
To characterize urban geometry, SPARTACUS-Urban takes building plan area fraction (λ p ) and the normalized building perimeter length (L) as a function of z for the city area of interest (hereafter "domain"). Thus, L(z) is the total building perimeter normalized by the horizontal area of the domain (units m −1 ). The SPARTACUS-Urban model discretizes the vertical profile into n layers where layer i has thickness z i and normalized perimeter length L i . Following this, the normalized wall area (A W ) is which is proportional to the frontal area index (λ f ) or projected wall area for a particular azimuthal direction (Raupach and Shaw 1982;Grimmond and Oke 1999;Sützl et al. 2020).
The SPARTACUS-Urban model assumes that walls face in all azimuthal directions with equal probability. Based on analysis of the geometry in real cities (Hogan 2019b), SPARTACUS-Urban assumes that the probability distribution of wall-to-wall horizontal separation distances, p ww (x), and the distribution of ground-to-wall separations, p gw (x), each follow an exponential distribution: where x is the horizontal wall-to-wall distance in any azimuthal direction, and X is the mean wall-to-wall distance, or 'e-folding' distance (Hogan 2019a, b). This exponential distribution allows the direct and diffuse streams of radiation to be attenuated according to the Beer-Lambert law (Hogan 2019a).
The assumption that buildings are randomly distributed horizontally has two other important consequences. First, the horizontal distribution of radiation between buildings (or vegetation, if included) need not be explicitly simulated; rather the horizontal-mean radiation field in each direction is computed as a function of height alone. Second, at a given height, the rate of interception of radiation by buildings is proportional the total building perimeter. However, in real cities Eq. 2 tends to underpredict the frequency of large building separations such as parks, parking-areas, and plazas (e.g., Fig. 6, Hogan (2019b)). Hence, the penetration of direct shortwave radiation to street level tends to be underpredicted when the sun is low in the sky, leading to an overprediction in absorption of shortwave radiation by walls. To address this, we replace Eq. 2 with the sum of two exponentials, with a weighting between them (G ww ): where X 1 and X 2 are the e-folding distance of each exponential. The weighting function (Eq. 3) better predicts the frequency of large building separations by up to three times that of Eq. 2, in theory improving the prediction of radiative fluxes. Applying Eq. 1 of Hogan (2019b) leads to an equation for p gw (x) with both the same form as Eq. 3 and the same X 1 and X 2 values, but a different weighting coefficient G gw To use the two-exponential model within SPARTACUS-Urban, we modify L(z) so that it varies with θ 0 , where θ k is the zenith angle of a stream k, where k = 0 indicates the solar zenith angle, defining an effective normalized perimeter length,L, given bŷ This is identical to Eq. 8 of Hogan (2019b) except for the use of an effective e-folding building separation (X ) in place of the e-folding building separation used to characterize the horizontal distribution of urban geometry. To derive an analytical relation betweenL and the two-exponential fit coefficients (G gw , X 1 and X 2 ), we assume that if all buildings had the same height, H, the fraction of direct solar radiation penetrating to street level (F 0g ) can be predicted exactly using Eq. 3 of Hogan (2019b) where x 0 = H tan θ 0 . As buildings are typically not all the same height, we use the mean building height (H ). Substituting our Eq. 3 (but for p gw ) in Eq. 6 gives However, as SPARTACUS-Urban follows an exponential building distribution (Eq. 2) we apply Eq. 6, leading to This is equal to Eq. 20 of Hogan (2019b), but using the effective e-folding building separation,X . Combining this with Eq. 5 giveŝ whereL 0 isL at the surface. This can be applied to any city with building footprint data, but probability distributions can only be computed using the Hogan (2019b) method near the surface for building densities, λ p > 0.01. UsingL 0 we scale L(z) at each height using the building cover fraction at that height This leads to the scaling factor to L(z) (Eq. 10) having a greater impact near the surface where λ p is larger, and a reduced impact as buildings thin out towards the top of any urban canopy, whereL tends toward L. The appropriateness of the single-and two-exponential methods in describing urban environments is discussed in Sect. 5, and their use in calculating fluxes is assessed in Sect. 6.3.

Modification to Treat Building Concavity
From SPARTACUS-Urban's use of L to describe the rate of interception of radiation by building walls, it follows that the width of the shadow cast by any building (W S ), averaged over all azimuthal illumination directions, is assumed to be equal to P/π where P is the building's perimeter length. This is true for convex shapes such as cylinders and cuboids, but not for many real-world buildings, which have fine structure in their perimeters, and so the shadows cast in SPARTACUS-Urban tend to be wider than reality. We therefore define the "concavity parameter" C as the ratio of the true perimeter length to the perimeter length of the equivalent convex hull, which for an individual building is given by Values of C are found to be 1 or greater, and to vary with height (Appendix 1, Fig. 12 and Online Resource 7). The effective concavity parameter of all the buildings at a particular height can be calculated by replacing the numerator and denominator of Eq. 11 each by the average over all buildings. In this work, the median of C at all heights above H for each domain is used as a scaling factor to L(z) at all heights (i.e., L(z)/C).

Description of the Discrete Anisotropic Radiative Transfer Model
The three-dimensional DART model (Gastellu-Etchegorry et al. 2015;Landier et al. 2018) simulates radiation propagation for heterogenous scenes that can include vegetation, buildings, a within-canopy atmosphere, and variations in ground height (i.e., topography). The latter can be imported using 3D vector models. The DART model has been evaluated using observations and other 3D radiative transfer models for vegetation (Sobrino et al. 2011;Widlowski et al. 2015), and has been applied in urban areas (e.g., Landier et al. 2018;Chrysoulakis et al. 2018;Morrison et al. 2020a). Radiative fluxes are calculated iteratively, with radiation tracked and emitted along a number of discrete directions within angular cones (Gastellu-Etchegorry et al. 2015) using a 3D array of voxels to facilitate radiation tracking. Radiation interacts with the scene elements in each voxel. Per-voxel scattered, absorbed, emitted, upwelling and downwelling radiative fluxes are updated after each iteration, with upwelling and downwelling fluxes for each voxel stored in the top face of each voxel. Here, we use DART's ability to calculate the radiative budget to assess the SPARTACUS-Urban emission and absorption of shortwave radiation.

Model Domains
Four types of urban form (F) ( Table 1) are used in the evaluation, from simplest to most complex: (1) Regular array of cubes that repeat on a regular grid (F REG1 and F REG2 , Table 1). Cubes are often used to approximate urban processes (e.g., Aida 1982;Kondo et al. 2005;Kanda et al. 2005;Kanda 2007; Morrison et al. 2018) as they create a regular grid street pattern that occurs in many cities in the U.S. and China (Figueiredo and Amorim 2007;Han et al. 2020).   Table 1 (2) Random cuboids (F RAND ) where building centroids are randomly located within a domain with random building heights, widths, and orientations. Twelve domains are used, but four are focussed on, with building fraction at the surface (λ p,0 ) and H values spanning those found in areas of real cities (F RAND1 -F RAND4 , Table 1) based on prior studies (Loridan and Grimmond 2012). This form type tests situations where the SPARTACUS-Urban building layout assumptions are met. These might be more typical of European cities, where street orientation is more random.
(3) Low level-of-detail (LOD) real-world geometry using building footprint data, with one height per building creating flat roofs and walls that do not taper, and flat ground ( Fig. 1a, d). Building footprints used are for part of central London (F Lon,L , Table 1). (4) High LOD real-world geometry where heights can vary across a building (Fig. 1b, e).
Parts of two cities are analyzed: a dense European megacity London (F Lon,H )-and an open low-density U.S. grid-city-Indianapolis (F Ind,H ).
The three real-world domains (i.e., 3 and 4) are 2000 × 2000 m 2 , to sample a wide range of streets with different widths, orientation, intersections, parking areas, plazas, and parks.
The DART model uses vector 3D building models to describe the urban form. For the low LOD, a raster digital surface model (DSM) and digital elevation model (DEM) are used to determine the building roof and ground level from the "Virtual London" building footprint dataset (Evans et al. 2006), using the 25th percentile of the DEM height, and the 75th percentile of the DSM height. Each building is assigned one height value from these building footprints. The 3D building models in the high LOD domains are created using the Morrison et al. (2020a, b) method from Google Earth imagery (Google Inc. 2019) and building footprints (Evans et al. 2006;Heris et al. 2020).
For SPARTACUS-Urban, profiles of λ p and L are calculated from rasterized (1 m resolution) building heights derived from the 3D building models. As SPARTACUS-Urban assumes no topographic variation within an individual NWP grid cell (i.e., flat), the DART 3D array of voxels is re-gridded to give heights relative to local ground level for high LOD scenes (Morrison and Benjamin 2021). To balance computational time and simulation resolution, the DART voxel resolution used is 1 m vertically and 15 m horizontally for both the realworld and F RAND scenes. For F REG , a vertical resolution of 0.5 m is used. In all scenes, SPARTACUS-Urban uses the same vertical resolution as DART.
As the real-world 3D building models are found to not conserve energy in DART, primarily because of periodic boundary conditions, the energy loss (always < 2%) needs to be redistributed. The rationale and method are explained in Appendix 2.
Fluxes from the DART model for F REG (Table 1) scenes are offset by the voxel vertical resolution, as DART provides the fluxes at the 'top face' of each voxel and all roofs are at 5 m, so at the top of a voxel. As SPARTACUS-Urban outputs wall absorption profiles between height intervals, DART wall absorption for F REG scenes is offset by 0.25 m. All data used and code are archived at https://zenodo.org/10.5281/zenodo.5145851.

Sun Angles and Albedos
Both the DART and SPARTACUS-Urban models require solar zenith angle, θ 0 , to be provided for a simulation. For computational simplicity we use three: directly overhead (0°, although unrealistic when out of the tropics), 45°, and low-sun conditions (75°). Incoming radiation at the top of the canopy is assumed to be directly from the sun; diffuse incoming radiation is set to zero. Similarly, a material albedo (α) is needed. We use two values: low (0.1) as observed in dense urban areas (e.g., 0.11, Kotthaus and Grimmond 2014) and high (0.5) as typical of 'cool' materials (e.g., Santamouris 2014; Santamouris et al. 2018;Jandaghian and Akbari 2018).
As SPARTACUS-Urban assumes that the azimuthal orientation of buildings is random, solar azimuth angle (Ω) is not specified by the model, whereas for DART the value of Ω is specified. Thus, DART has varying shadow patterns with Ω used. The DART simulations use four values for the simpler F RAND and F REG cases, and eight (at 45°intervals) for the real-world cases. The final DART fluxes for comparison use the mean across all Ω intervals.

Evaluation Statistics
To quantify SPARTACUS-Urban performance, we compare SPARTACUS-Urban and DART profiles of: mean shortwave upwelling (SW ↑ ) and downwelling clear-air (SW ↓ ) fluxes, and mean wall (a Wall ) and roof shortwave absorption (a Roof ). The fluxes have units of watts per square metre (W m −2 ) of the entire horizontal scene (rather than per square metre of the clear-air region excluding buildings), while the absorptions have units of W m −3 , since we divide the absorption in a layer by the layer thickness to obtain a resolution-independent quantity. Thus, the vertical integral of a Wall and a Roof provide the total wall and roof absorptions (again per unit area of the entire horizontal scene). Unlike the vertical walls assumed by SPARTACUS-Urban for all domains, for the high LOD geometry in DART there is no simple way to distinguish or define roofs and walls. Hence, we combine the wall and roof absorption (a Wall+Roof ) for the evaluation of the high LOD scenes. All fluxes and absorptions are normalized by the bottom of atmosphere (BOA) shortwave flux (SW ↓,BOA ). This is defined as the incoming shortwave flux across a horizontal plane above the tallest roughness elements in a scene (Gastellu-Etchegorry et al. 2015;Wang et al. 2020).
Profiles of a Wall and a Roof are compared using the normalized mean-absolute error (nMAE) and normalized mean-bias error (nMBE) expressed as a percentage of the mean DART where a SU and a DART are the flux at each height from SPARTACUS-Urban or DART, respectively. The metrics nMAE and nMBE are calculated at scene resolution (e.g., 1 m) vertically from 1 m to the maximum height in DART (H max ). The SW ↑ flux profiles are evaluated using the normalized bias error at a specified height (nBE), expressed as a percentage of the DART flux: The scene albedo is evaluated using SW ↑ at the top of the canopy (H max in DART). We also use the metric nBE to evaluate the total ground absorption (a Ground ).

Test of the SPARTACUS-Urban Geometry Assumptions
We examine the underpinning SPARTACUS-Urban assumption-that urban buildings are randomly distributed, or equivalently their horizontal separations follow an exponential distribution (Eq. 2)-by analyzing probability distributions from real cities and domains containing randomly placed cuboids (Table 1).
For the high density F RAND3 domain, the 'true' probability density of wall-to-wall (p ww ) and ground-to-wall (p gw ) separations (Fig. 2) are calculated following Hogan (2019b) with 1 × 1 m 2 resolution building rasters, analyzed in four azimuthal directions 45°apart. Both the p ww and p gw distributions fit a single-exponential well (Fig. 2b, c) for separations up to 200 m, indicating F RAND3 satisfies SPARTACUS-Urban's assumption of randomly distributed buildings. This behaviour is seen for all F RAND domains. Here we use Eq. 2, where X is obtained from Eq. 5 with the surface value of L (denoted L 0 ).  Table 1) within a 2 × 2 km 2 domain with a plan area fraction at the surface, (λ p,0 ) of 0.5 and a mean height (H ) of 7 m: a Plan view, randomly placed cuboid buildings and probability density with a single-exponential (Eq. 2) fits to the b wall-to-wall and c ground-towall probability distributions Fig. 3 For a 2 × 2 km 2 area of central London at low LOD (F Lon,L ): a wall-to-wall and b ground-to-wall probability distribution, with single-(Eq. 2, blue) and two-exponential fit (Eq. 3 where X 1 = 21.6 m, X 2 = 57.4 m and G gw = 0.534, red), and c corresponding effective normalized perimeter length at the surface as a function of solar zenith angle,L 0 (θ 0 ), as predicted by Eq. 9 when applied to the actual data (black) and fitted (red) ground-to-wall probability distributions, and the true perimeter length at the surface (blue) here eight azimuthal directions are used to determine p gw and p ww , but offset by 22.5°to not align with major streets orientations (e.g., north-south or east-west in Indianapolis, Fig. 1c). Comparison of the calculated probability distribution to both the single-and twoexponential fits for central London indicates that the latter is a better fit for F Lon,L (Fig. 3) and F Lon,H (Online Resource 1). For F Lon,L , the single-exponential fit diverges from the p ww distribution at approximately 200 m and from the p gw distribution at approximately 100 m (Fig. 3b, c), whereas for F Ind,H (Fig. 4) the single-exponential fits diverge at slightly greater distances (these numbers increasing to around 300 and 200 m, respectively). By contrast, the two-exponential (Eq. 3) predicts the larger building-separations much better than the single-exponential in both cities. For F Lon,L the effective normalized building perimeter length, L 0 , decreases when θ 0 > 30°(black line Fig. 3c). This is computed using Eq. 5 with the true ground-to-wall probability distribution (i.e., Fig. 3b), the mean building height (H = 25.5 m), and F 0g (Eq. 6). This shows that more direct solar radiation reaches the surface when the sun is low in the sky (i.e., less chance of wall interception) compared to purely randomly distributed buildings. Using the actual normalized perimeter of 0.055 m −1 (blue, Fig. 3c), equivalent to the single-exponential assumption, would be expected to lead to an overpredicted interception of direct solar radiation by walls for larger θ 0 . The L value obtained from Eq. 9 with the two-exponential method agrees well with L obtained using the true probability (Fig. 3c), and the same behaviour is seen for Indianapolis (Fig. 4c) (H = 17.9 m). Thus, we expect that the two-exponential fit should improve SPARTACUS-Urban simulations for real-world cities. This is tested in Sect. 6.

Regular Cubes
Comparison of shortwave radiative flux profiles simulated with DART and SPARTACUS-Urban (single-exponential, Eq. 2) in a low-density regular array of cubes (F REG1 ) shows that SW ↓ decreases closer to the surface when the zenith angle θ 0 = 75° (Fig. 5) because more radiation is intercepted by buildings. Hence, less shortwave radiation penetrates to ground level. For all θ 0 values, the roof absorption, a Roof , remains constant, with nMAE = 0 (i.e., machine precision). As the buildings are all the same height, a Roof can be computed exactly from the building fraction and albedo. Maximum values of a Wall increase with θ 0 (Fig. 5c, g, k), with a Wall increasing with height due to building shadowing at the surface (as θ 0 increases). Azimuth angle (Ω) variations change shadow patterns and alter the wall area exposed to shortwave radiation. Hence, the a Wall vertical profiles differ between DART and SPARTACUS-Urban. The DART model's fluxes are averaged across four Ω values (Sect. 4.2). The SPARTACUS-Urban model's a Wall profiles are within the DART range arising from Ω (Fig. 5c, g, k shading) when θ 0 = 0°and 45°, but not when θ 0 = 75°and α = 0.1 near the surface (approximately 1 m). Errors in a Wall are lowest when θ 0 = 45°(nMAE between 9.9 and 17%, nMBE between −7.5 and 1.8%). For the scene albedo the nBE are < 1.2%, with values highest if the sun is overhead (θ 0 = 0°). When α = 0.1, nBE in SW ↑ and SW ↓ increase, but are still < 2% (Online Resource 2). These results are better than expected, given the grid arrangement of the buildings do not have an exponential distribution of building separations (i.e. as SPARTACUS-Urban assumes, Sect. 2.1).
The larger plan area index of F REG2 (Table 1) causes SW ↓ to decrease more as height decreases (for high θ 0 ) (Fig. 6a, e, i). The form F REG2 also increases mutual building shadowing, reducing the shortwave radiation penetrating to the surface. With less shortwave radiation escaping, the scene albedo decreases with increasing θ 0 . The metric nBE is larger (up to 10%, Table 2b) cf. F REG1 . Maximum values of a Wall increase as θ 0 increases (Fig. 6c, g, k). However, a Wall decreases more rapidly as height decreases than in F REG1 , due to the increased shadowing from the buildings/cubes. The value of nMAE is larger (cf. F REG1 ) for a Wall (up to 35%, Table 2b). The peak in a Roof remains at 5 m, as all buildings are of equal height. The SPARTACUS-Urban model generally overpredicts the SW ↑ and SW ↓ profiles,  Table 2 Evaluation of SPARTACUS-Urban (relative to DART) for F REG scenes (Table 1) (a) F REG1 and (b) F REG2 , for three solar zenith angles (θ 0 : 0, 45, 75°) and one albedo (α: 0.5). Upwelling and downwelling clear air shortwave flux profiles (SW ↑ and SW ↓ ) assessed with the normalized bias error (nBE, Eq. 14) and wall and roof absorption (a Wall , a Roof ) profiles assessed with the normalized mean absolute error (nMAE, Eq. 12) and the normalized mean-bias error (nMBE, Eq. 13) and underestimates a Wall at the top of the canopy (Fig. 6) in F REG2 . Similarly to F REG1 , when α = 0.1 nMAE in a Wall can be up to 35% when θ 0 > 45°(Online Resource 2).

Random Cubes
Four F RAND domains (F RAND1 to F RAND4 ) are intended to test SPARTACUS-Urban performance across the λ p,0 and H extreme combinations, with more results for eight other F RAND domains given in Online Resource 3. All F RAND simulations use the single-exponential model (Eq. 2), as it fits the building distribution data well (Fig. 2). Figure 7 shows the agreement between DART and SPARTACUS-Urban for each θ 0 and α for profiles of a Wall , a Roof , and SW ↓ for F RAND3 . Overall, for F RAND1 -F RAND4 , the nBE and nMAE are less than 6% for all quantities (Table 3), as SPARTACUS-Urban's urban form assumptions (Sect. 5) are fulfilled. The SPARTACUS-Urban model agrees better with DART when λ p,0 and H are small, as buildings are further apart so there is less within-canyon scattering and building shadowing. The largest differences between DART and SPARTACUS-Urban are seen for a Wall between 1 and 5 m for θ 0 = 45°, 75°. The SPARTACUS-Urban model underestimates SW ↑ for θ 0 = 0°and 45°. In F RAND1-2 , when λ p,0 = 0.05, nBE < 0.7% compared to nBE = 3.4-5.0% when λ p,0 = 0.5 (F RAND3-4 ). When λ p,0 = 0.05, nMAE < 1%, except for a Wall in F RAND1 , where nMAE = 2.1%. Although performance becomes poorer as λ p,0 increases, F RAND3-4 errors do not exceed 5.5% when θ 0 = 75°and α = 0.5. The differences in nBE and nMAE magnitudes are larger with an increase in λ p,0 (F RAND1 -F RAND3 ), compared with an increase in H (F RAND1 -F RAND2 ). This is also seen in the additional scenes in Online Resource 3.

Real-World Geometry
For the real-world urban form in SPARTACUS-Urban, both the single-(Eq. 2) and twoexponential (Eq. 3) fits are used, allowing assessment of the impact of the building layout Building height distribution profiles have spikier patterns for F Lon,L than F Lon,H (blue and orange, Online Resource 4c), with the largest differences between 25 and 50 m. These spikes occur because individual buildings in the F Lon,L domain each have only one height, and are aggregated per 1 m interval. Despite this, the vertical profiles of a Roof in SPARTACUS-Urban and DART are still close (Fig. 8d, h, l).
The SPARTACUS-Urban model has good agreement to DART for F Lon,L vertical flux profiles (Fig. 8), although the agreement is generally poorer with increasing θ 0 . The SPARTACUS-Urban model always underestimates SW ↑ , with nBE values within 7% of DART (Table 4a). In the SPARTACUS-Urban model, a Ground is overestimated but with nBE < 6%. The SPARTACUS-Urban model is generally in better agreement to the DART model for both SW ↑ and a Ground when using the two-exponential method . This is most evident as θ 0 increases. Neither nMAE nor nMBE exceed 7.3% for a Wall . Generally, the SPARTACUS-Urban model underestimates a Roof , with nMAE < 12%, and nMBE up to −4.3% (increases to 13% and −6.3% for the single-exponential method). When α = 0.1 for the two-exponential, the maximum magnitudes of nBE for SW ↑ increases (10%, Online Resource 5) but nBE for a Ground , and nMAE for a Roof are similar (cf. α = 0.5).
The vertical absorption profiles (a Wall , a Roof , a Wall+Roof ) for the London scenes are well captured by SPARTACUS-Urban (Figs. 8 and 9). The a Wall+Roof maxima between DART and SPARTACUS-Urban (Fig. 9f, i) disagree mainly because of the need to adjust for intra-scene local topography heights in DART, in contrast to SPARTACUS-Urban where the ground is assumed to be flat. The SPARTACUS-Urban model fluxes for F Lon,H are within 12% of DART, Table 3 As Table 2, but for F RAND (Table 1) Table 4b). Both scene albedo and transmission to the surface (Table 4b) are underestimated at all θ 0 . SPARTACUS-Urban overestimates the a Wall+Roof profiles (Fig. 9c, f, i) leading to less shortwave radiation at ground level (with reduced SW ↓ with decreasing height as θ 0 increases). This leads to less shortwave reflected to the top of the canopy, reducing the scene albedo. The largest differences between the single-and two-exponential results occurs when θ 0 = 75° (Table  4b). Simulations for F Lon,H using α = 0.1 (cf. 0.5) have higher nBE magnitudes for both SW ↑ (3.0-5.7%), and a Ground (7.5-36%, Online Resource 5). Error metrics (nBE, nMAE, nMBE) are notably larger when using the single-exponential.

Comparison to the Single-Layer Infinite-Street-Canyon Assumption
Given the current urban models within NWP models commonly assume an urban form consisting of an infinite canyon with buildings of the same height and flat roofs, we assess SPARTACUS-Urban relative to one model of this type, Harman et al. (2004). This solves a Table 4 As Table 2 Harman et al. (2004) longwave radiation by modifying the configuration, so assumptions are met for both, viz.: buildings all the same height, and exponential model of urban geometry (Hogan 2019b). Hogan (2019a) found excellent agreement between SPARTACUS-Urban using eight streams, supporting the use of the discrete ordinate method for urban radiative transfer.
Here, for shortwave radiation, we compare SPARTACUS-Urban and Harman et al. (2004) against DART for cases when the assumptions of the two models are not necessarily satisfied. The Harman method is used with its usual configuration (i.e., exchange coefficients consistent in a single-layer infinite street canyon as in Sect. 3 of Hogan (2019b)), and we implement the 2 × 2 matrix inversion approach of Harman et al. (2004), as outlined in Sect. 4.2 and Eq. 4 of Hogan (2019a). This approach assumes two parallel infinite length buildings have constant height, H, separated by street of constant width, W , with a fixed H/W .
Care is taken to ensure that in all comparisons, the total area of ground, wall, and roof is equal between the three models. For the Harman simulations, the height H is set equal to H 320 M. A. Stretton et al.  (Table 1), and the building fraction equal to λ p,0 . The value of H/W is calculated using the operational method in Eq. 3 of Hertwig et al. (2020), where λ f is calculated for each domain using Eq. 1 with the true wall area of the domain calculated from z and L(z). From Eq. 15 we obtain W using H . For SPARTACUS-Urban, we use the two-exponential form. Analysis is undertaken for both random cuboid and real-world scenes.
Unlike SPARTACUS-Urban, the Harman et al. method cannot predict vertical profiles, so the comparison of wall and roof absorptions is limited to vertically integrated quantities. Values of SW ↑ at the top of the canopy are calculated for the H max in each scene for DART and SPARTACUS-Urban. These are compared using nBE (Eq. 14). We expand on results for an albedo of 0.5 here, with results of the comparison for an albedo of 0.1 in Online Resource 6. Run times for five SPARTACUS-Urban configurations are compared. These have varying numbers of diffuse streams per hemisphere (N = 1 to 8) and layers: n = 1 (i.e., as Harman et al. (2004)) to 6 (e.g., reasonable for operational NWP), to 151 (i.e., this work). The computer time of a single-threaded run for a SPARTACUS-Urban profile with the most basic configuration (N = 1, n = 1) is fast (12 µs) but six times longer than for the Harman model (Table 5). The F Lon,L scenes (Sect. 6.3; SPARTACUS-Urban configuration: N = 8, n Table 5 Absolute run-time of the three models (Harman, SPARTACUS-Urban, and DART) for the low LOD London domain (F Lon,L , Table 1)  = 151) have a much longer run time (9.2 ms) but SPARTACUS-Urban is~2.5 million times faster than DART despite DART using 14 parallel threads (Table 5).
Overall, Harman et al. (2004) has the best agreement with DART SW ↑ (nBE < 6.3%, Fig. 1) but generally overestimates a Wall , a Roof , and a Wall+Roof . The best agreement between Harman et al. (2004) and DART is found for F RAND scenes (Fig. 11), however, SPARTACUS-Urban performs better (cf. Harman). Harman absorption errors increase with θ 0 , with nBE up Fig. 11 Comparison of SPARTACUS-Urban, Harman et al. (2004) and DART values for real-world scenes at low and high LOD (F Lon, L , F Lon,H , F Ind,H ), and random cuboid scenes (F RAND ) for three solar zenith angles (θ 0 : 0°, 45°, 75°) and an albedo of 0.5: upwelling clear air flux at the top of the canopy (SW ↑ ), and total wall, roof, and ground absorption (a Wall , a Roof , a Ground ). For high LOD scenes a Wall and a Roof are combined (a Wall+Roof ) for evaluation. Numbers on each bar are the nBE (Eq. 14) between the SPARTACUS-Urban model/Harman approach and the DART model. Error bars for DART span the range between if no correction for the energy imbalance is made, and if double this correction is made (Appendix 2) to 18% or 32% (a Wall and a Roof respectively), compared with nBE for SPARTACUS-Urban (up to 13.7% or 0.8%).
For real-world domains, SPARTACUS-Urban performs better than the Harman approach, particularly for low θ 0 (Fig. 11). The nBE SPARTACUS values are in the range 1.3-14.3% for a Roof , a Wall , and a Wall+Roof , while nBE Harman is 0.3-31% for the same quantities. The value of a Roof is overestimated by both models in all cases, except for the SPARTACUS-Urban model for F Lon,L at 45°and 75°(magnitudes of nBE SPARTACUS < 4.3%, Fig. 11, cf. nBE Harman < 21.8%), which is expected for the Harman approach as roof shadowing is neglected. Both models predict SW ↑ well (nBE all < 7%), but the SPARTACUS-Urban model almost always performs better (Fig. 11). The worst performance for both models is for a Ground , with both models generally underestimating it (nBE of 0.4-26.6% for SPARTACUS-Urban and 1.2-56.4% for the Harman method). Overall, using SPARTACUS-Urban has a smaller nBE (cf. Harman) when evaluated using DART. For scenes where nBE Harman (cf. DART) are lower than nBE SPARTACUS (i.e., better performance), the differences in nBE are < 5%.

Conclusions
Evaluation of the multi-layer SPARTACUS-Urban shortwave fluxes is undertaken using reference calculations from the explicit 3D radiative transfer model DART. The SPARTACUS-Urban model computes the vertical profiles of fluxes and absorption rates in urban scenes, which is crucial for vertically resolved urban energy balance models. A range of urban geometries were considered: regular arrays of cubes, cuboids with random placement and heights, to real city complexity (London and Indianapolis).
The SPARTACUS-Urban approach performs well when the SPARTACUS-Urban assumption of randomly distributed buildings is fulfilled. This is particularly evident for low building densities (λ p,0 = 0.05) where the normalized bias error (nBE) and normalized mean absolute error (nMAE) < 1%. The SPARTACUS-Urban model and the DART model agree less well as building fraction increases (λ p,0 = 0.5, nBE and nMAE < 5.5%). The largest nBE and nMAE occur when the solar zenith angle is highest (θ 0 = 75°). For all random cuboid scenes presented, all nBE and nMAE are below 6%.
The shortwave radiative fluxes for real cities (London and Indianapolis) have nBE magnitudes of less than 7% for effective scene albedo, and nBE generally less than 15% for ground absorption. Exceptions to the latter occur for the high level of detail London domain when θ 0 = 45°and 75°. Errors (nMAE) for the wall and roof absorption (low LOD) are less than 7% and 12%, respectively. The combined wall and roof absorption (high LOD domains) is always overestimated by SPARTACUS-Urban (nMAE < 15%), which leads to underestimation in the effective albedo of the scene, and underestimation in the transmission of shortwave radiation at the surface. However, the structure of the vertical profiles of fluxes and absorptions are captured well by SPARTACUS-Urban. Overall, upwelling profiles are best predicted by SPARTACUS-Urban. For the low LOD London domain, shortwave downwelling is typically overestimated, and roof absorptions are generally underestimated, in contrast to the high LOD domains. The performance of SPARTACUS-Urban in Indianapolis is slightly worse than for London scenes, which could be related to the grid-like street layout being further from the random building distribution assumed by SPARTACUS-Urban. Nonetheless, the Indianapolis domain used here still contains parks and diagonally oriented streets, which make the domain and building separations sufficiently random enough that SPARTACUS-Urban still performs well.
Regular cube arrays tested have a similar form to earlier urban radiation studies. The roof absorptions modelled in SPARTACUS-Urban are exact (cf. DART), as all buildings have the same height. The smallest differences between SPARTACUS-Urban and DART are found in scene albedo and ground absorption, where nBE < 2.2%. This increases to < 18% for denser cube arrays. Across all scenes, the largest differences between DART and SPARTACUS-Urban are found in wall absorption, with nMBE between 1.8 and −22% (nMAE 9.9% and 31%) depending on cube density. These errors in wall absorption are greater than in the real scenes, which given the regular spacing does not meet the randomly distributed buildings SPARTACUS-Urban assumption, this result is not unexpected. It is plausible the low nMAE and nBE results may be associated with the open cube spacing reduces building shadowing effects.
A modification to the original SPARTACUS-Urban method is introduced here, relaxing the strict assumption that the distribution of building separations fits a single-exponential (Eq. 2). This is replaced with a two-exponential method, allowing for effective building edge length to vary with θ 0 (Eq. 3). This new method is proposed because Eq. 2 underpredicts the frequency of large building separations, and thus underpredicts the fraction of solar radiation reaching ground level. Using the two-exponential method both improves the predicted probability distributions (cf. 'true' distributions for London and Indianapolis) and reduces the SPARTACUS-Urban model's radiative flux errors by up to a factor of a half. A further correction is applied to both the single-and two-exponential methods to account for the concavity of real buildings, leading to better representation of the width of building shadows. The range of possible concavity values in real-world cities means that SPARTACUS-Urban simulations could be calibrated to fit DART simulations. For NWP, the concavity value uncertainty for an individual domain is smaller than the uncertainty in L(z) itself.
There is scope to further refine SPARTACUS-Urban, including adjusting building shadowing. As SPARTACUS-Urban assumes the shadow cast by a building falls onto neighbouring buildings, or on gaps between buildings, is proportional to the roof area, this means shadows are randomly overlapped with roofs that they fall on. However, buildings often have roofs at low and high heights that are effectively next to each other when viewed from overhead. So, higher parts of a building can shadow lower roofs, rather than the street-level. Correcting this could improve the SPARTACUS-Urban performance further.
Comparison to the Harman et al. (2004) single-layer infinite street canyon model for randomly distributed cuboid, and real-world geometries found it performs best (c.f. DART) for random cuboid scenes. However, SPARTACUS-Urban generally performs better and notably in real-world scenes. The model results are most similar in their effective scene albedo predictions (nBE generally < 7%). The Harman approach overestimates the roof absorption, which is expected as the single-layer infinite street approach neglects roof shadowing.
Overall, our results show the SPARTACUS multi-layer approach to modelling radiative transfer in urban areas agrees well with the more complex and computationally demanding radiative transfer model DART when modelling real-world cities. The SPARTACUS-Urban model is the first multi-layer urban radiation model to achieve this, whilst being computationally cheap enough to be incorporated into weather and climate models (Hogan 2019a). This work surpasses previous evaluations that compare radiative transfer models to small scale observations or to more simplistic radiative transfer models (e.g., Harman et al. 2004;Krayenhoff and Voogt 2007;Aoyagi and Takahashi 2012;Krayenhoff et al. 2014). The single-exponential distribution that underpins SPARTACUS-Urban performs well but can be improved by using a two-exponential method. It is not yet certain if the extra complexity of implementing this two-exponential is justified, given the uncertainty in urban morphology datasets. Such datasets are required to compute required model inputs for SPARTACUS-Urban to describe the urban form (i.e., vertical descriptions of the urban canopy), which would be required if SPARTACUS-Urban is to be applied into a large-scale model.
Further investigation is needed to ascertain the amount of data required to describe building geometry worldwide, and how this impacts radiative fluxes. Further evaluation should be completed with SPARTACUS-Urban in the longwave, as this is a significant term in the urban surface energy balance (Oke 1988). As SPARTACUS-Urban can be integrated within existing urban surface energy budget models, the results of this shortwave evaluation provide a promising start to improving the treatment of the complex urban structure in NWP models.

Supplementary Information
The online version contains supplementary material available at https://doi.org/ 10.1007/s10546-022-00706-9. For the London domains, C ranges between 1.04 and 1.66 at each height interval, with smaller values for the low (F Lon,L ) than high (F Lon,H ) LOD domain (Fig. 12). Generally, C decreases with height in real cities (Fig. 12, Online Resource 7), because individual buildings in an area of the city may be short and wide, or tall and narrow, and so cast different sized shadows. Cross-sections of the high LOD London domain support this, with low height levels, where the building density is high, having a larger C (cf. higher heights where smaller, taller building occur) (Online Resource 7d). Given this, we calculate C at all heights, and use median C from all heights above H , applying this as the scaling factor to L(z) (black, Fig. 12). Values below H are excluded given both the presence of courtyards, and that taller buildings will shadow shorter buildings (rather than vice versa).
We examine the impact on the shortwave fluxes, with C from 1.1 to 1.7 but constant with height (Sect. 4.3). Across the three solar zenith angles (Fig. 13), the C value has the greatest impact on the SW ↑ profiles for the F Lon,H domain, and less impact on the SW ↓ and a Wall+Roof profiles. Fig. 13 As Fig. 9, but SPARTACUS-Urban profiles indicating the variation from using a concavity parameter of 1.1 (solid), 1.4 (dotted), and 1.7 (dashed). SPARTACUS-Urban uses the two-exponential (Eq. 3) only

Appendix 2: Redistribution of Lost Energy in the DART Model Output
We calculate energy imbalance (E, W m −2 /W m −2 ) for each DART run, using where the total wall, ground and roof absorption are given by a Wall , a Ground and a Roof , respectively and SW ↑,top is the shortwave upwelling clear air flux at the top of the canopy. E = 0 for a perfectly conserving model. The DART energy imbalance is always less than 2%, and usually less than 1% (Fig. 14b). Generally, energy conservation is highest for overhead sun conditions (θ 0 = 0°) and decreases with geometry complexity to be lowest for the real-world domains (F Lon , F Ind ). F REG (Table  1) domains have a lower voxel resolution, so do not follow this relation, as DART energy conservation tends to improve as voxel resolution increases (e.g. 1 m → 0.5 m).
Processes identified that contribute to energy loss are (Fig. 15): (1) rays passing 'underground' of the domain when there are holes in the 3D building model, rays hitting internal walls of buildings Both can occur from the domain periodic boundary conditions, if buildings and topography are cut at the domain edge (Wang et al. 2020). The F RAND , F Lon,H and F Ind,H ( Table  Fig. 14 Energy loss (Eq. 16) in DART scenes (Table 1) for three solar zenith angles (θ 0 : 0°, 45°, 75°) for albedos: a 0.1, and b 0.5, for all F REG , F RAND1-4 , and real-world (F Lon,L , F Lon,H , F Ind,H ) scenes tested (Sect. 4.1) 1) 3D building models contain some buildings that cross the domain edge. As the domain cuts though it creates 'open' buildings allowing rays to pass directly underneath the building footprint model or interact with the 'internal' walls of the buildings (Fig. 15b). With non-flat topography (e.g., F Lon,H and F Ind,H ) the periodic boundary conditions create gaps through which rays can pass (Fig. 15a). F REG scenes have neither of these features, whereas F RAND scenes have only building split, hence energy loss increases with increasing domain complexity. As the real-world scenes have the largest energy loss, the DART missing energy is assumed to be lost through the building walls and the ground, attributed to processes (1) and (2) above.
For high LOD scenes, we redistribute lost energy into the walls (Fig. 15b, d) first as a ratio of total A W to total a Wall . The approximated (A W,edges ) wall area missing at the edge of the domain is used to calculate the average wall absorption a Wall,edges = A W , edges a Wall We note A W,edges is an overestimate, as building walls in the repeated units may overlap. We assume a Wall,edges is the amount of energy absorbed by these extra walls. The amount of total energy available for exchange by the wall processes (E Wall ) in Fig. 15d is equal to a Wall,edges /(1 − α). This energy is split between SW ↑ , a Wall , and a Ground (Fig. 15d). The remainder of the lost energy is attributed to the ground process (E Ground ) (Fig. 15a), and is redistributed (Fig. 15c) into a Ground , a Wall , and SW ↑ , where the two-thirds of the radiation Fig. 15 Energy lost in DART occurs at the edge of periodic model domains (grey) if there is a unmatched topography, and/or b buildings missing external walls. To distribute the energy through the c ground (E Ground ), and d walls (E Wall ), where the purple arrows (i) represent the extra energy into the ground and walls, and the red (ii) and green (iii) represent the first and second reflections of shortwave radiation, to determine the total extra shortwave radiation into the walls, ground, and upwelling (a Wall,extra , a Ground,extra , and SW ↑,extra ) reflected from the ground is distributed into a Wall . Combining these processes (Fig. 15) leads to the total added energy.
This is distributed at all height levels using a scaling factor (e.g., a Wall /a Wall,extra ). For F Lon,L , the buildings at the edge of the domain have complete walls with flat topography, so we distribute lost energy through both the wall and ground processes (Fig. 15) equally, assuming E Wall = E Ground . Although we redistribute the energy through these processes, we note this may not match the true DART results if external walls and topography are corrected. The uncertainty in these numbers is less than the range if solar azimuth angle is varied, and so it not shown in the vertical flux profiles compared in Sect. 6. In Sect. 7, this uncertainty is shown by an error bar across: the fluxes if no correction for the energy imbalance is made, and the fluxes if double the correction is made.