Benchmarking a new urban scheme in the ORCHIDEE 2.2 land surface model

Urban areas change natural surface energy and water balances, yet some land surface models still represent cities as natural surfaces. In this study, we present the development of a one-tile urban scheme for the ORCHIDEE land surface model, designed to improve the representation of urban processes, particularly for high-resolution applications. The scheme incorporates key urban parameters such as albedo, building height, thermal properties, and imperviousness. We propose a novel physically-based approach for representing imperviousness, by modifying saturated hydraulic conductivity to account for both surface and subsurface impacts. Off-line simulations across 20 urban flux tower sites show improved performance in sensible and latent heat fluxes with the new urban scheme, compared to the original baresoil representation. Mean absolute error (MAE) evaluation confirms improved model skill, aligning with benchmark results from the Urban-PLUMBER intercomparison. The scheme also captures expected urban hydrological signatures, such as increased surface runoff, though some reductions in subsurface runoff may be less consistent with observed urban recharge patterns. This work lays the foundation for applying ORCHIDEE at the basin scale and in high-resolution convection-permitting simulations for urban hydroclimate studies. Perspectives include refining thermal parameter choices, integrating anthropogenic heat fluxes, and conducting high-resolution simulations to assess hydrological performances using observed streamflow data.

1 Introduction

Land Surface Models (LSMs) simulate surface processes such as the surface energy and water budgets and are designed to be coupled to atmospheric models, with an emphasis on robust land–atmosphere exchanges across diverse environments (Blyth et al., 2021). Historically, their development followed two pathways shaped by the application scale. GCM-oriented LSMs, such as ORCHIDEE, CLM, and JSBACH, operate on coarse grids (tens to hundreds of kilometers; e.g. Fisher and Koven, 2020) and focus on global water and carbon cycles (Krinner et al., 2005; Lawrence et al., 2019; Raddatz et al., 2007). By contrast, LSMs used in mesoscale and limited-area models, such as Noah-MP, TERRA, and ISBA-SURFEX, have long used finer resolutions, favoring more accurate landscape representation and parameterizations relevant at kilometer scale over restricted domains or shorter periods (Doms et al., 2011; Masson et al., 2013; Niu et al., 2011), but neglecting groundwater processes (Barlage et al., 2021). Increasing computing power is now blurring this distinction, as GCMs move toward kilometer-scale grid spacing (Prein et al., 2026; Segura et al., 2025; Stevens et al., 2019), and mesoscale models increasingly run in climate-mode over large domains and long periods (Lucas-Picher et al., 2021; Schär et al., 2020). This convergence highlights the need for scalable, efficient parameterizations valid across a wide range of environments. Urban land cover exemplifies this challenge: while explicitly represented for decades in mesoscale models (e.g., Masson, 2000), it remains relatively new in GCM-oriented LSMs where coarse grids previously limited its relevance. The shift toward kilometer-scale modeling and rapid urban expansion (Güneralp et al., 2020; Seto et al., 2012), motivate the development of physically based urban schemes in these models. Yet, urban processes are still insufficiently represented in GCM-oriented LSMs despite their strong influence on land–atmosphere (Langendijk et al., 2024; Michau et al., 2025) and land–subsurface interactions (Bhaskar et al., 2016; Oudin et al., 2018).

To reproduce the characteristic signatures of the urban surface energy budget, an urban representation must account for several key surface properties and fluxes that differ markedly from natural areas: (1) radiative properties, particularly albedo and longwave-radiation trapping; (2) aerodynamic properties, including surface roughness and urban geometry; (3) thermal properties, especially heat capacity and thermal conductivity; (4) vegetation cover, imperviousness, and water availability; and (5) anthropogenic heat emissions. First, cities modify the absorbed shortwave radiation through changes in the effective urban canopy albedo relative to surrounding rural areas because of the materials used in cities and the geometry of buildings and streets (Kotopouleas et al., 2021; Qin, 2015). Urban geometry also affects longwave exchanges by trapping infrared radiation within the canopy. Together with the high heat capacity and thermal conductivity of urban materials, these radiative effects contribute to increased heat storage and often higher urban surface temperatures. The sensible heat is modified by cities buildings, which alter wind speed (depending on the street geometries can either increase or decrease), and increase surface roughness (Sun et al., 2024). The sensible heat is also modified by the different surface temperatures (usually higher than surrounding natural areas) of urban materials which change the temperature gradient between the surface and the air (Barlow, 2014). By reducing the vegetation cover, and by increasing impervious surfaces, cities typically decrease the latent heat flux (Loridan and Grimmond, 2012). This reduction in latent heat leaves more available energy at the surface to be dissipated as sensible heat. In urban environments, ground heat flux has been estimated to account for 30 %–40 % of net radiation, in contrast to just 10 % in rural areas (Rigo and Parlow, 2007). This larger heat storage results from both the higher heat capacity and thermal conductivity of urban materials and the greater surface area available for heat exchange and storage within the urban canopy (Bernabé et al., 2015; Johra, 2021; Wonorahardjo et al., 2020). This refers to what is typically called the storage heat in urban climate studies, with urban materials storing significant heat during the day and releasing it at night (Lindberg et al., 2020). Even if not explored within this manuscript, the anthropogenic heat flux (QF) is a term typically added to urban surface energy budgets. It refers to anthropogenic heat emissions from buildings, from domestic and non-domestic sources (human activities such as heating, cooling, transportation, and industrial processes), with diurnal, seasonal, and regional patterns (Hamilton et al., 2009).

Urban representations in LSMs began to evolve in the early 2000s, moving from simple to increasingly process-based descriptions. Early approaches treating urban surfaces as natural land covers failed to capture distinct urban radiative, momentum, and heat fluxes at higher spatial resolutions (Lipson et al., 2024). The first step toward explicitly representing cities in LSMs was the introduction of urban bulk schemes (Myrup, 1969; Oke, 1988), in which grid cell properties such as albedo, roughness length, and surface thermal characteristics were adjusted to reflect urban conditions while still treating the surface as a single homogeneous entity. To improve the representation of the urban energy balance, and to address street scale climate questions, canyon-based models were implemented in LSMs. The Town Energy Budget (TEB; Masson, 2000) implemented an explicit urban canyon geometry within a LSM, distinguishing roofs, walls, and roads and allocating separate energy budgets to each surface. This was followed by multi-layer configurations, such as the Building Effect Parametrization (BEP; Martilli et al., 2002), which added vertical discretization of urban elements and enabled more refined exchanges of heat, momentum, and moisture. The most recent classification of geometries of urban schemes (Lipson et al., 2024) places them on a continuum: (i) non-urban representations (bare soil, vegetation), (ii) one-tile “slab” schemes that aggregate urban surfaces, (iii) two-tile schemes distinguishing between roofs and roads, (iv) canyon-based schemes resolving roofs, walls, and roads separately usually with multi-layers in the canyon, and (v) more complex three-dimensional approaches.

Beyond morphology, urban schemes in LSMs diverge in the processes and fluxes they represent. Comparisons have shown that latent heat flux is often the least well represented, and historically even neglected, under the assumption that vegetation and water availability are minimal in cities (Best and Grimmond, 2015). Over the last decade, however, evidence has emphasized the importance of vegetation for correctly simulating latent heat fluxes and model performance during wet conditions. Current urban schemes differ in how vegetation is represented: it may be excluded entirely, handled through a sub-tile vegetation fraction with separate energy and water budgets, or integrated directly into the canyon geometry (Krayenhoff et al., 2014; Lemonsu et al., 2012; Redon et al., 2020).

On top of vegetation, hydrology has also been shown as a key element to capture correctly latent heat fluxes. Urbanization modifies surface and subsurface hydrological pathways by increasing soil sealing, compaction, changing preferential flow pathways, and through water use infrastructures (Bhaskar et al., 2016; Fletcher et al., 2013). These changes generally increase surface runoff and shorten hydrological response times, while their effects on soil moisture, groundwater recharge, low flows, and evapotranspiration depend on local soil properties, vegetation, water use, drainage connectivity, and sewer or stormwater networks (Bhaskar et al., 2016; Fletcher et al., 2013; Oudin et al., 2018; Saadi et al., 2020). This complexity is partly reflected in the use of the term “imperviousness”, which is not unique in urban hydrology. Depending on the modeling and spatial context, it can refer to the total impervious surface fraction, the fraction effectively connected to the drainage network, or an effective hydrological control on infiltration and surface runoff generation (Jacobson, 2011; Saadi et al., 2020). This distinction is important for ULSMs because urban hydrology influences not only surface runoff, but also soil water availability and evapotranspiration, which directly links the water and energy budgets.

In ULSMs, urban hydrology has been introduced progressively through several modelling choices. Early urban canopy schemes primarily focused on radiative, thermal, and aerodynamic exchanges, while hydrological processes were absent or simplified, with artificial surfaces commonly represented as impervious and water exchange mainly limited to surface interception, evaporation from wet sealed surfaces, and runoff (Kusaka et al., 2001; Martilli et al., 2002; Masson, 2000). Later developments introduced more explicit hydrological distinctions between urban surface types. One direction was to distinguish hydrological behaviour between pervious and impervious urban fractions, for example through pervious and impervious canyon-floor fractions in CLM-U (Oleson et al., 2008). TEB-related developments add infiltration through artificial surfaces, and connected rainfall interception on roofs and roads toward a soil component (Lemonsu et al., 2007). Other schemes refined the representation of water stored on impervious surfaces and its evaporation after rainfall, as in TERRA-URB (Wouters et al., 2015). Hydro-microclimate approaches further extended urban hydrology to include subsoil beneath built surfaces, lateral water transfer between pervious and impervious tiles, sewer drainage, and irrigation, as in TEB-Hydro (Stavropulos-Laffaille et al., 2021, 2018). Despite these advances, recent intercomparison work shows that urban water-balance representation remains uncertain in ULSMs, especially regarding water-balance closure, storage dynamics, and runoff parameterization (Jongen et al., 2024). Progress is also constrained by observations, because urban sites providing turbulent fluxes do not provide all water-balance terms needed to evaluate runoff, irrigation, drainage, and storage changes consistently at the same spatial scale (Jongen et al., 2024).

The Urban-PLUMBER intercomparison project (Lipson et al., 2024) demonstrated that simpler models in terms of morphology, when evaluated as a group, often outperformed more complex ones in simulating local-scale fluxes for large-scale applications such as regional or global climate modeling. The results from the Urban-PLUMBER intercomparison also indicate that model performance is not strictly linked to geometric complexity: schemes that represent hydrology, vegetation, and imperviousness realistically, regardless of whether they are slab, two-tile, or canyon, can perform well across turbulent heat fluxes, whereas complex schemes often underperform when these elements are simplified (Lipson et al., 2024). These findings underscore that accurately representing water availability, imperviousness, and vegetation may be as important as the level of geometric detail.

In its current configuration, the ORCHIDEE (Organising Carbon and Hydrology In Dynamic Ecosystems) land surface model represents urban areas as natural surfaces, more specifically, as bare soil. This simplification limits the model's ability to capture key urban processes, including sensible and latent heat flux partitioning and urban runoff generation. As ORCHIDEE is increasingly applied at fine spatial resolutions, including in convection-permitting coupled simulations below 3 km (Shahi et al., 2022), the absence of an explicit urban representation becomes a critical limitation. ORCHIDEE already includes a detailed and physically based hydrological scheme and a sub-tile vegetation representation through plant functional types (PFTs). However, because urban areas are currently treated as bare soil, ORCHIDEE lacks a mechanism to represent impervious surfaces and to assign surface parameters that reflect urban form and materials.

To address the limited representation of urban processes in ORCHIDEE, particularly the effects on energy and water surface budgets, we propose the development of a one-tile urban scheme within the model. This approach builds on advancements in urban parameterizations already implemented in other land surface models. Our goal is to design a scheme that can readily incorporate urban characteristics datasets while remaining computationally efficient. Our objectives are twofold. First, we design an urban one-tile scheme that captures the dominant controls on the urban surface energy and water budgets while remaining consistent with ORCHIDEE's existing structure and computational constraints. We parameterize key urban characteristics (e.g., albedo, thermal conductivity, heat capacity, building height) such that they can be assigned from global or local urban datasets. Second, we introduce a parsimonious, physically motivated, representation of imperviousness. Rather than resolving separate water budgets for pervious and impervious urban sub-tiles, the proposed approach preserves the existing ORCHIDEE soil hydrology with an urban soil tile where imperviousness is represented through its effect on the hydraulic properties of the soil column. Specifically, imperviousness modifies saturated hydraulic conductivity as a function of depth, allowing it to influence infiltration, surface runoff generation, soil water storage, and vertical water transfer. This formulation is designed for one-tile, kilometer-scale and coarser LSM applications, where the objective is to represent their aggregate effect on the water budget and surface fluxes feedback to the atmosphere. We evaluate the new scheme at urban and suburban flux tower sites, compare its performance to existing ORCHIDEE configurations, and assess the relative importance of imperviousness-related hydrological processes for urban energy and water fluxes. Section 2 describes the model and the new developments. Section 3 presents the experimental setup and evaluation strategy. Section 4 provides the results, Sect. 5 discusses their implications, and Sect. 6 summarizes the conclusions and perspectives.

2 Model Development

2.1 ORganizing Carbon and Hydrology In Dynamic EcosystEms (ORCHIDEE)

ORCHIDEE is the Land Surface Model developed at IPSL (Institut Pierre-Simon Laplace) and used in global and regional climate models, as well as in standalone mode. It simulates coupled water, energy, and carbon fluxes between the land surface and the atmosphere (Boucher et al., 2020; Krinner et al., 2005). This study uses ORCHIDEE version 2.2, revision 8133; further details and access information are provided in the Code and data availability section.

As shown in Fig. 1, land cover is represented within each grid cell as a mosaic of plant functional types (PFTs). The standard configuration includes up to 15 PFTs: one bare-soil PFT, eight high-vegetation PFTs (forest types), and six low-vegetation PFTs (including grasslands and croplands). Each PFT is associated with its own radiative, morphological, physiological, and phenological characteristics. Water and energy fluxes associated with vegetation and natural soils including transpiration, evaporation of water intercepted by the vegetation canopy, and soil evaporation, are computed at the PFT level and are then aggregated to the grid cell scale using a flux-aggregation approach (calculated separately for each PFT and then aggregated to the grid-cell scale as area-weighted fluxes). In contrast, energy fluxes associated with radiative and aerodynamic parameters (e.g., albedo, roughness length) are directly calculated at grid-cell scale using PFT-weighted parameters.

Figure 1 Schematic representation of the energy- and water-budget calculations in ORCHIDEE for a single grid cell, its plant functional types (PFTs), and the three associated soil columns representing bare soil, high vegetation, and low vegetation. Yellow arrows denote fluxes calculated directly at the grid-cell scale, green arrows denote fluxes calculated separately for each PFT and subsequently aggregated to the grid-cell scale, and dark-blue arrows denote fluxes calculated independently for each soil column. For example, sensible heat flux (H) is calculated at the grid-cell scale, evapotranspiration (ET) at the PFT scale, and surface runoff and subsurface runoff at the soil-column scale. Each soil column is discretized into 12 hydrological layers extending to a depth of 2 m, while the thermal discretization extends to 18 m. Abbreviations: P, precipitation; ET, evapotranspiration; H, sensible heat flux; SW, shortwave radiation; LW, longwave radiation; PFT, plant functional type.

Each grid cell contains three soil columns, one for bare soil, one for high vegetation, and one for low vegetation. These columns have the same depth, vertical discretization, soil texture, and texture-dependent hydraulic properties. They differ in their vegetation-related processes: the high- and low-vegetation columns have distinct rooting profiles and root-water uptake, whereas the bare-soil column has no root extraction. Their soil-moisture states and water budgets therefore evolve independently. Each soil column is vertically discretized (12 layers, down to a depth of 2 m for hydrology, and 18 m for soil temperature computation), allowing the model to simulate temperature and moisture transfer through the soil profile. For each soil column, ORCHIDEE calculates the surface runoff, and subsurface runoff, representing gravitational outflow at the bottom of the soil column. The organization of the PFT mosaic and associated soil columns, together with the spatial scale at which the different fluxes are calculated, is summarized in Fig. 1.

2.1.1 Representation of the energy fluxes and surface energy budget

The surface energy budget is often described as:

$$\begin{array}{}\text{(1)}& R_{\mathrm{net}}=H+\mathit{LE}+G\end{array}$$

where R_net is the net radiation, representing the total balance between incoming shortwave (SW) and longwave (LW) radiation and their respective outgoing components, H is the sensible heat flux, LE is the latent heat flux, and G is the ground heat flux. The heat stored and released by the soil or urban substrate is not neglected, but is represented through (G), which is used as the upper boundary condition of the one-dimensional heat diffusion equation for the soil temperature profile. Thus, no additional storage term is introduced in the surface energy budget.

ORCHIDEE computes these terms following its internal energy budget formulation. Building upon previous PhD and accreditation to supervise research (HDR) theses (Arboleda-Obando, 2023; Arjdal, 2023; Chéruy, 2018), we summarize here only the components relevant to the development of an urban one-tile scheme. For clarity, snow-covered conditions and floodplain processes, although represented in ORCHIDEE, are omitted from the description below because they are not relevant to the urban one-tile developments presented here. At each time step, at the grid cell scale, ORCHIDEE computes a surface temperature that is consistent with the radiative, turbulent, and conductive fluxes. This surface temperature is then used to compute the surface energy fluxes at the current time step and is carried forward, together with the updated soil temperature profile, as part of the thermal state used at the following time step.

The net radiation term is computed as:

$$\begin{array}{}\text{(2)}& R_{\mathrm{net}}=\left(\mathrm{1}-\mathit{\alpha }\right){\mathrm{SW}}_{\mathrm{in}}+\mathit{ϵ}{\mathrm{LW}}_{\mathrm{in}}-\mathit{ϵ}\mathit{\sigma }{T_{\mathrm{surf}}}^{\mathrm{4}}\end{array}$$

where α is the albedo (unitless); σ is the Stefan–Boltzmann constant ($\mathrm{5.67}\times {\mathrm{10}}^{-\mathrm{8}}$ W m⁻² K⁻⁴); ϵ (unitless) is the surface longwave emissivity, fixed to 1 in these simulations, corresponding to a blackbody approximation in the thermal infrared; and T_surf is the surface temperature (K). Equation (2) therefore defines the radiative term of Eq. (1). It is not independent from the other fluxes because T_surf also controls the turbulent and conductive fluxes.

Each vegetation PFT has its own value of albedo (α_leaf), and the value of albedo for baresoil (α_MODIS) is read from MODIS (Moderate Resolution Imaging Spectroradiometer) data (Lurton et al., 2020). The use of MODIS-derived bare-soil albedo fields accounts for the spatial variability of soil reflectance resulting from heterogeneities in soil composition, moisture, and texture. Unlike vegetation, where albedo can be more homogeneously modeled based on physiological properties, bare soil requires an empirical approach. MODIS satellite data provides precise measurements of these variations, ensuring that local differences in soil reflectance are accurately captured. The resulting albedo of a grid cell is a mean of albedos of bare soil and of the PFTs constituting the grid cell, weighted by the fractional areas of the PFTs of the grid cell:

$$\begin{array}{}\text{(3)}& \mathit{\alpha }=f_{{\mathrm{PFT}}_{\mathrm{1}}}{\mathit{\alpha }}_{\mathrm{MODIS}}+\sum _{\mathrm{2}\text{–}\mathrm{15}}^i{f}_{{\mathrm{PFT}}_i}{\mathit{\alpha }}_{{\mathrm{leaf}}_i}\end{array}$$

where $f_{{\mathrm{PFT}}_{\mathrm{1}}}$ is the baresoil fraction of the grid cell; $f_{{\mathrm{PFT}}_i}$ is the PFT_i (vegetation PFTs) fraction of the grid cell; ${\mathit{\alpha }}_{{\mathrm{leaf}}_i}$ is the albedo of each vegetated PFT.

For each grid cell at each time step, ORCHIDEE solves the energy budget equations based on a bulk-diffusion approach (Barella-Ortiz et al., 2013). It allows the calculation of the sensible and latent heat fluxes exchanged between the land surface and the atmosphere accounting for the various components of the latent heat flux: soil evaporation, plant transpiration, evaporation of water intercepted by the vegetation canopy (interception loss), snow sublimation. The sensible and latent heat fluxes are calculated respectively thanks to Eqs. (4) and (5).

$$\begin{array}{}\text{(4)}& {\displaystyle }H=\mathit{\rho }v{C}_{\mathrm{d}}{c}_p\left(T_{\mathrm{surf}}-T_{\mathrm{air}}\right)\text{(5)}& {\displaystyle }\mathit{LE}=L_{\mathrm{evap}}\mathit{\beta }\mathit{\rho }v{C}_{\mathrm{d}}\left(q_{\mathrm{surf}}-q_{\mathrm{air}}\right)\end{array}$$

where ρ is air density (kg m⁻³); v is the horizontal wind speed (m s⁻¹); C_d is the drag coefficient (unitless); c_p is the specific heat capacity of dry air (1004.675 J kg⁻¹ K⁻¹); L_evap is the latent heat of evaporation or sublimation (2.5008×10⁶ or 2.8345×10⁶ J kg⁻¹); β is the limiting factor of potential total evapotranspiration (unitless); T_surf and T_air are the surface and air temperature respectively (K); and q_surf and q_air are the specific moisture at the surface and in the air (kg kg⁻¹), respectively. Equations (4) and (5) define the H and LE terms of Eq. (1). They are linked to Eq. (2) through T_surf as the same surface temperature controls outgoing longwave radiation, the surface–air temperature gradient, and the saturation humidity entering the latent heat flux calculation.

Turbulent transfer of the fluxes between the surface and the atmosphere, is represented thanks to the drag coefficient C_d. For each PFT in the standalone version of ORCHIDEE, C_d which can also be expressed as $\mathrm{1}/(r_{\mathrm{ah}}\cdot u)$, where u is the wind speed and r_ah the aerodynamic resistance, is calculated as follows:

$$\begin{array}{}\text{(6)}& C_{\mathrm{d}}={\displaystyle \frac{k^{\mathrm{2}}}{\ln {\left(\frac{z_{\mathrm{lev}}}{h_{\mathrm{PFT}}\cdot z_{\mathrm{0}}}\right)}^{\mathrm{2}}}}\end{array}$$

where k is the Von Karman constant; z_lev is the height of the first atmospheric layer (10 m in standalone simulations); h_PFT is the height of the canopy of the considered PFT (fixed) and z₀ is the dimensionless roughness-length-to-height ratio and is fixed to 0.0625 ($\mathrm{1}/\mathrm{16}$, default value in ORCHIDEE for all PFTs in standalone simulations). Thus, the product h_PFT⋅z₀ in Eq. (6) corresponds to the aerodynamic roughness length. The aerodynamic roughness length is therefore not prescribed as an additional independent parameter but is diagnosed from the PFT height. For the urban PFT, h_PFT corresponds to the prescribed site-specific building height; consequently, its aerodynamic roughness varies among sites. Because ORCHIDEE calculates a single energy budget per grid cell, the PFT-level drag coefficients are subsequently aggregated to the grid-cell scale using the corresponding land-cover fractions.

The remaining term in Eq. (1) is the ground heat flux G. It is computed from the thermal gradient between the surface and the upper soil or substrate layer:

$$\begin{array}{}\text{(7)}& G={\displaystyle \frac{\mathit{\lambda }}{C_p}}{\displaystyle \frac{\partial T_{\mathrm{surf}}}{\partial z}}\end{array}$$

where λ (W m⁻¹ K⁻¹) and C_p (J kg⁻¹ K⁻¹) are respectively the soil thermal conductivity and the soil heat capacity.

This expression constitutes the upper boundary condition of the one-dimensional heat diffusion equation used to compute the evolution of soil temperature along the vertical profile. The lower boundary condition is a zero-heat-flux condition at the bottom of the soil column (18 m), following the standard ORCHIDEE configuration. Soil temperature is prognostically calculated by solving the Fourier heat diffusion equation:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{\mathit{\lambda }}{C_p}}{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial z^{\mathrm{2}}}}\end{array}$$

where T is the soil temperature (K).

Finally, the soil thermal properties also change according to the soil water content of each soil layer:

$$\begin{array}{}\text{(9)}& {\displaystyle }{C}_p=\left(C_{\mathrm{wet}}+\mathrm{1}-C_{\mathrm{dry}}\right)S+C_{\mathrm{dry}}\text{(10)}& {\displaystyle }\mathit{\lambda }=\left({\mathit{\lambda }}_{\mathrm{s}}\cdot {\mathit{\lambda }}_{\mathrm{w}}-{\mathit{\lambda }}_{\mathrm{dry}}\right)S+{\mathit{\lambda }}_{\mathrm{dry}}\end{array}$$

where, λ_s and λ_w are the heat conductivities for solid soil and water, respectively (W m⁻¹ K⁻¹); λ_dry is the heat conductivity of dry soil (W m⁻¹ K⁻¹); S is the total soil saturation degree (m³ m⁻³); C_dry and C_wet are respectively the dry and saturated soil heat capacity (kJ m⁻³ K⁻¹). The hydrological and thermal discretizations are identical in the upper 2 m of the soil column. Within this common domain, hydrological soil-moisture variables are transferred to the thermal grid, including interpolation of volumetric water content to the thermal-layer interfaces. Below 2 m, soil moisture is not prognostically simulated. Instead, the volumetric water content is assumed to be vertically uniform and equal to that of the deepest hydrological layer. The total water content assigned to each thermal layer therefore scales with its thickness. This assumption is consistent with the gravitational-drainage boundary condition at the bottom of the hydrological column.

2.1.2 Representation of water fluxes at the grid cell scale

The description of the computation of the water budget and its different fluxes is detailed in the technical documentation “The hydrol module of ORCHIDEE” (Ducharne et al., 2018), as well as in PhD theses using the ORCHIDEE model (Campoy, 2013; d'Orgeval, 2006). The following sections take up the description of this report and these PhD theses, focusing on the relevant information for the development of an urban one-tile scheme.

The surface water budget reflects the balance between water inputs and outputs at the surface. In ORCHIDEE, it is defined by the partitioning of precipitation into evapotranspiration, surface and subsurface runoff, and a change in soil moisture. Each grid cell contains three soil columns for water budget calculations: one for bare soil, one for high vegetation, and one for low vegetation. These columns are 2 m deep and discretized into 12 layers whose thickness increases with depth. The vertical redistribution of soil moisture, here volumetric water content θ (m³ m⁻³), in these columns follows a diffusion-based scheme that relies on the Richards equation (De Rosnay et al., 2002), with surface and subsurface runoff acting as the upper and lower boundary conditions, respectively and is stated as:

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial \mathit{\theta }\left(z,t\right)}{\partial t}}=-{\displaystyle \frac{\partial J\left(z,t\right)}{\partial z}}-s\left(z,t\right)\end{array}$$

where z is depth below the soil surface (m), t is time (s), J is the flux field (m s⁻¹), and s is the sink term (m³ m⁻³ s⁻¹). The flux J follows Darcy's law (Darcy, 1856), later extended to unsaturated soils by Buckingham (1907).

$$\begin{array}{}\text{(12)}& J\left(z,t\right)=-D\left(\mathit{\theta }\left(z,t\right)\right){\displaystyle \frac{\partial \mathit{\theta }\left(z,t\right)}{\partial z}}+K\left(\mathit{\theta }\left(z,t\right)\right)\end{array}$$

where K(θ) and D(θ) are the hydraulic conductivity and diffusivity (in m s⁻¹ and m² s⁻¹ respectively). The relationships between K, D, and θ are given by the Mualem (1976)–Van Genuchten (1980) model:

$$\begin{array}{}\text{(13)}& {\displaystyle }K\left(\mathit{\theta }\right)=K_{\mathrm{s}}\sqrt{{\mathit{\theta }}_{\mathrm{f}}}{\left(\mathrm{1}-{\left(\mathrm{1}-{\mathit{\theta }}_{\mathrm{f}}^{\mathrm{1}/m}\right)}^m\right)}^{\mathrm{2}}\text{(14)}& {\displaystyle }D\left(\mathit{\theta }\right)={\displaystyle \frac{\left(\mathrm{1}-m\right)K\left(\mathit{\theta }\right)}{m\mathit{\mu }}}{\displaystyle \frac{\mathrm{1}}{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}}{\mathit{\theta }}_{\mathrm{f}}^{-\mathrm{1}/m}{\left({\mathit{\theta }}_{\mathrm{f}}^{\mathrm{1}/m}-\mathrm{1}\right)}^{-m}\end{array}$$

where K_s is the saturated hydraulic conductivity (m s⁻¹), μ (m⁻¹) corresponds to the inverse of the air entry suction, m is a dimensionless parameter, related to the classical Van Genuchten parameter n, by $m=\mathrm{1}-\mathrm{1}/n$; θ_r is the residual water content (m³ m⁻³), and θ_f the relative humidity (m³ m⁻³) defined as ${\mathit{\theta }}_{\mathrm{f}}=(\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}})/({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}})$, with θ_s the saturated water content (m³ m⁻³).

In version 2.2 of ORCHIDEE, the vertical profile of saturated hydraulic conductivity K_s(z) is prescribed as a function of depth. Its value at any depth z results from three components:

• i.

  Reference saturated hydraulic conductivity ($K_{\mathrm{s}}^{\mathrm{ref}}$) at the depth z_lim = 0.30 m, derived from the dominant USDA (U.S. Department of Agriculture) soil texture of the grid cell (12-class texture classification). Each grid cell's dominant soil texture influences soil transfer (Ducharne et al., 2018; Tafasca et al., 2020).

• ii.

  An exponential decrease with depth, controlled by the dimensionless factor F_K(z). This decrease applies from z_lim to the bottom of the soil column (2 m).

• iii.

  A root-induced enhancement in the upper soil layers, represented by a factor F_Kroots(z,c), where c denotes the soil-column index. This enhancement is applied only to vegetated soil columns, whereas F_Kroots(z,c) is set to 1 for the bare soil column, which therefore retains a constant K_s over the first 30 cm (d'Orgeval et al., 2008; De Rosnay et al., 2002).

The resulting expression for the saturated hydraulic conductivity profile is:

$$\begin{array}{}\text{(15)}& K_{\mathrm{s}}\left(z\right)=K_{\mathrm{s}}^{\mathrm{ref}}{F}_K\left(z\right)F_{\mathrm{Kroots}}\left(z,c\right)\end{array}$$

where F_K(z) controls the exponential decrease with depth:

$$\begin{array}{}\text{(16)}& F_K\left(z\right)=\min \left(\max \left(\exp \left(-f\left(z-z_{\lim }\right),{\displaystyle \frac{\mathrm{1}}{F_K^{\max }}}\right),\mathrm{1}\right)\right)\end{array}$$

Apart from the texture-dependent $K_{\mathrm{s}}^{\mathrm{ref}}$, the other three parameters (z_lim, f, $F_K^{\max }$) are constants, i.e. they do not depend on the PFT, the soiltile, and the soil texture. The depth z_lim is set to 0.30 m, the decay factor f is set to 2 m⁻¹, and the dimensionless parameter $F_K^{\max }$ is set to 10. Therefore, F_K(z) cannot be smaller than $\mathrm{1}/F_K^{\max }=\mathrm{0.1}$.

Infiltration is computed using a Green–Ampt type approach (Green and Ampt, 1911). When precipitation exceeds the infiltration capacity at the surface, the non-infiltrated fraction becomes surface runoff (Horton or infiltration-excess runoff). Surface runoff is computed independently for each soil column, and the grid cell surface runoff is obtained by summing the contributions from the three columns. At the bottom of the 2 m soil column, ORCHIDEE applies a free gravitational drainage boundary condition by default. This means that water leaves the soil whenever hydraulic conductivity at the bottom layer is nonzero. This downward flux is referred to as subsurface runoff. As for surface runoff, grid cell subsurface runoff is the sum of the subsurface runoff from the three soil columns. The total runoff is the sum of subsurface runoff and surface runoff. There are no subsurface lateral fluxes between adjacent soil grid cells in the model.

In the ORCHIDEE model, evapotranspiration includes evaporation from soil, and vegetation canopy, along with plant transpiration. Evapotranspiration (E_total) is the link between the energy balance and the water balance, and can be written as:

$$\begin{array}{}\text{(17)}& E_{\mathrm{total}}={\mathit{\beta }}_{\mathrm{ET}}{E}_{\mathrm{pot}}\end{array}$$

and:

$$\begin{array}{}\text{(18)}& E_{\mathrm{pot}}=\mathit{\rho }v{C}_{\mathrm{d}}\left(q_{\mathrm{surf}}-q_{\mathrm{air}}\right)\end{array}$$

The actual evapotranspiration E_total (kg m⁻² s⁻¹) is thus calculated from a potential evapotranspiration E_pot (kg m⁻² s⁻¹), and a stress function β_ET (unitless). This stress function β_ET has its own value for each type of evaporation:

$$\begin{array}{}\text{(19)}& {\mathit{\beta }}_{\mathrm{ET}}={\displaystyle \frac{A_{\mathrm{v}}\left({\mathit{\beta }}_{\mathrm{1}}+{\mathit{\beta }}_{\mathrm{2}}\right)+A_{\mathrm{g}}{\mathit{\beta }}_{\mathrm{3}}}{A_{\mathrm{L}}}}\end{array}$$

where β₁, β₂ and β₃ are the individual stress functions on interception loss, transpiration and bare soil evaporation respectively, and A_v, A_g and A_L are the fractions of vegetated fraction, bare soil fraction, and total land area, respectively. β₁ = 1, since the local evaporation flux from the intercepted water proceeds at its potential rate. The local stress functions β₂ and β₃ are ≤1, and become closer to zero as soil moisture gets limiting.

Once the latent heat flux is calculated with Eqs. (5) and (17), ORCHIDEE calculates the different terms contributing to total evapotranspiration:

$$\begin{array}{}\text{(20)}& E_{\mathrm{total}}=E_{\mathrm{inter}}+E_{\mathrm{trans}}+E_{\mathrm{soil}}\end{array}$$

with each term depending on the different β previously mentioned, with for instance the E_soil being calculated as:

$$\begin{array}{}\text{(21)}& E_{\mathrm{soil}}={\mathit{\beta }}_{\mathrm{3}}\mathit{\rho }\left|\mathit{V}\right|C_{\mathrm{d}}\left(q_{\mathrm{surf}}-q_{\mathrm{air}}\right)\end{array}$$

with L_evap the latent heat of vaporization.

2.2 New urban parametrization for ORCHIDEE v2.2

2.2.1 Adding an Urban PFT with prescribed albedo and height

The first step in implementing the urban one-tile scheme was to introduce an Urban PFT into ORCHIDEE's land-cover representation. Urban-specific parameter values are prescribed for this PFT. For the energy budget, we assign urban values of albedo, and height of roughness elements. For the water budget, the key input parameter is the impervious surface fraction, ranging from 0 to 1. In this study, this fraction is not used to define separate pervious and impervious hydrological sub-tiles, nor does it represent the fraction directly connected to the drainage network. Instead, it is converted into an effective hydrological parameter that modifies saturated hydraulic conductivity as a function of depth, thereby affecting infiltration, runoff generation, soil-water storage, and vertical water transfer within the ORCHIDEE soil columns. All other PFT attributes remain identical to those of the bare-soil PFT, consistent with the slab-scheme assumption of no vegetation within the urban tile.

In ORCHIDEE, albedo is first computed at the PFT level and then aggregated to the grid cell scale by weighting each PFT contribution by its fractional cover (Eq. 3). For the bare-soil PFT, albedo is prescribed from MODIS, whereas vegetation PFTs use prescribed leaf albedo values ranging from 0.18 to 0.34, depending on vegetation type. For the new Urban PFT, albedo is handled in a flexible manner: it can either be prescribed as a single constant value specified by the user in the simulation namelist or read from an external urban-albedo dataset that ORCHIDEE interpolates onto the model grid. The same structure applies to the parameter controlling effective urban height, which influences aerodynamic roughness and turbulent fluxes. Urban height can correspond to a uniform value defined by the user or be derived from a gridded dataset describing building height or urban morphology. This approach ensures that the urban representation remains compatible with large-scale climate simulations while allowing the use of more detailed datasets when available.

2.2.2 Adding an urban soil column with imperviousness-dependent hydraulic conductivity

In a second step, we introduce a dedicated soil column associated with the Urban PFT, enabling soil water and heat processes to be represented separately from those of bare soil, high vegetation, and low vegetation. Water transfers into the soil depend on several factors, including the saturated hydraulic conductivity K_s (Eqs. 13 and 15). To represent the effects of imperviousness, K_s is reduced to limit water transfer into the soil column. To represent imperviousness within the urban soil column, we modify Eq. (15) by replacing the vegetation root factor F_Kroots(z,c) with a constant urban factor F_urb, which does not depend on depth:

$$\begin{array}{}\text{(22)}& K_{\mathrm{s}}^{\mathrm{urb}}\left(z\right)=K_{\mathrm{s}}^{\mathrm{ref}}{F}_K\left(z\right)F_{\mathrm{urb}}\end{array}$$

When F_urb equals 1, the original profile of K_s for bare soil remains unmodified. The lowest value assigned to F_urb is 0.1, as a value of 0 would imply no water transfer within the column, and values lower than 0.1 would result in unrealistically low saturated hydraulic conductivities for certain soil textures (typically with clay). To define the new factor F_urb, two options are considered. Option “Urban1” considers F_urb as a constant (same over all urban soil column), while in option “Urban2” F_urb depends on the imperviousness. The Urban1 option states F_urb equal to 0.1 to have $K_{\mathrm{s}}^{\mathrm{ref}}$ reduced to 10 % of its original value:

$$\begin{array}{}\text{(23)}& F_{\mathrm{urb}\mathrm{1}}=\mathrm{0.1}\end{array}$$

While the Urban2 option calculates F_urb as follows:

$$\begin{array}{}\text{(24)}& F_{\mathrm{urb}\mathrm{2}}=\mathrm{1}-\mathrm{0.9}{f}_{\mathrm{imp}}\end{array}$$

where f_imp is the value of imperviousness of the urban column and ranges from 0 to 1, so K_s(z) decreases proportionally to f_imp. When site-specific information is unavailable, the modeler may prescribe a representative value. For example, the Copernicus Urban Atlas distinguishes discontinuous medium-density urban fabric ($\mathrm{0.30}<f_{\mathrm{imp}}\le \mathrm{0.50}$), discontinuous dense urban fabric ($\mathrm{0.50}<f_{\mathrm{imp}}\le \mathrm{0.80}$), and continuous urban fabric (f_imp>0.80; European Environment Agency, 2020).

These two options are illustrated in Fig. 2, which represents the saturated hydraulic conductivity for a loam soil over the 2 m soil column. The orange line presents the original values from the baresoil column. The purple line presents the saturated hydraulic conductivity for the Urban1 hypothesis. The Urban1 hypothesis is designed to represent a very impervious conceptualization of cities. The pink line presents the saturated hydraulic conductivity for the Urban2 hypothesis, considering an imperviousness equal to 55 %.

Figure 2 Example profiles of saturated hydraulic conductivity for a reference loam soil (in orange) modified by the two options of imperviousness for an imperviousness equal to 0.55, and a reference saturated hydraulic conductivity of 249.6 mm d⁻¹, values as in Ducharne et al. (2018).

2.2.3 Assigning urban thermal properties to urban-dominated grid cells

The thermal conductivity (λ) and heat capacity (C_p) of the soil are calculated based on soil textures and soil water content (Eqs. 9 and 10), and are defined at the grid cell scale; they do not depend on PFTs. To take into account the specific urban thermal conductivity and heat capacity, we change the values of λ and C_p for the grid cells where urban areas cover 50 % or more of the grid cell. In grid cells where the urban fraction is below 50 %, the standard ORCHIDEE values calculated from soil texture and soil moisture are retained. This straightforward binary approach (where the grid cell is assigned either fixed urban values for thermal properties, or retains its initial values) was chosen because these variables are defined at the grid cell scale, rather than at the PFT or soil tile scale. The new prescribed values for urban areas are taken from the Noah-MP land surface model (He et al., 2023): thermal conductivity is set to 3.24 W m⁻¹ K⁻¹ and the heat capacity to 1890 kJ m⁻³ K⁻¹. These fixed values are not adjusted according to urban typology, building morphology, or surface-area density. In comparison, the range of values for grid cells with less than 50 % of urban cover is between 0.79 and 2.02 W m⁻¹ K⁻¹ with a median value of 1.20 W m⁻¹ K⁻¹ for thermal conductivity, and between 1191 and 2093 kJ m⁻³ K⁻¹ with a median value of 1632 kJ m⁻³ K⁻¹ for heat capacity. Urban areas therefore exhibit thermal conductivities that exceed those of most other soil textures.

3 Experiments and Methodology

3.1 Experiments

To evaluate the performance of the ORCHIDEE land surface model with the new urban parameterization, including urban albedo, building height, thermal conductivity, heat capacity, and imperviousness, a comprehensive series of 1D standalone simulations is conducted. These simulations are designed to systematically assess the model's sensitivity to these urban parameters in varying climatic and land-use contexts. The simulations are performed across 20 urban flux tower sites (Table 1), sourced from the harmonized gap-filled dataset of Lipson et al. (2022b), part of the Urban-PLUMBER project. For each site, the dataset provides all atmospheric forcing variables needed to force a one-dimensional land surface model, as well as site-level information on urban parameters and land use cover. These sites span 13 different countries and encompass a diverse range of climate zones: 11 sites are in temperate climates, 7 in continental climates, 1 in a tropical climate, and 1 in a dry climate. The Supplement provides a world map of these flux tower sites, including their climates and observed periods (Lipson et al., 2022b). Beyond climatic diversity, these urban sites also represent a wide spectrum of land-use characteristics, ranging from highly urbanized areas to low-density urban zones with significant vegetation cover (e.g., trees, grasslands, or croplands). Table 1 also provides a detailed breakdown of site-specific land-use characteristics. The Minneapolis station features two distinct datasets due to differing dominant wind influences, each capturing different land-cover compositions.

Table 1 Sites description with their sites land use and land cover characteristics and soil texture (Lipson et al., 2022b). The surfaces of the Urban, Tree, Low Veget, Bare soil, and Water land covers are expressed in land surface fraction. Clay and sand are expressed as soil-texture fractions, with the remaining fraction corresponding to silt. Building height is expressed in meters, and latitude (Lat) and longitude (Lon) are expressed in decimal degrees (°).

The numerical framework consists of multiple simulations for each urban site, each serving a specific purpose in the model evaluation process. The simulations are categorized into benchmark runs and urban parameterization experiments, as outlined in Table 2. The benchmark simulations are Baresoil and Lowveget. In all simulations, the observed non-urban land-cover fractions of each site, including trees, low vegetation, bare soil, and water, are retained. Only the fraction classified as urban is represented differently among the experiments. In Baresoil, the urban fraction is assigned bare-soil properties, corresponding to the default treatment of urban areas in ORCHIDEE. In Lowveget, the urban fraction is assigned the properties of the dominant low-vegetation type in the surrounding area, either grassland or cropland, providing a non-urban baseline. Three distinct urban configurations are also tested, all incorporating the newly developed urban albedo, building height, thermal conductivity, and heat capacity parameters. The key difference among them lies in the treatment of imperviousness: Urban0 which overlooks the influence of imperviousness on K_s to isolate the influence of energy-related parameters; Urban1 which incorporates the first imperviousness option (Urban1); and Urban2 which implements the second option (Urban2) for comparative assessment. Further details on grid cell albedo and building height used across all simulations and sites are provided in Table S1 in the Supplement. These values correspond to the parameter-aggregated approach described for albedo in Eq. (3). For each site, the model is first spun up for 40 years, using repeated cycles of 10 years of meteorological forcing provided by Urban-PLUMBER, until soil moisture and temperature reach equilibrium. The subsequent simulation period covers 12–16 years depending on the station (excluding spin-up). Fluxes are computed and saved on a 30 min time-step.

Table 2 Simulations description.

3.2 Evaluation strategy

To evaluate the impact of the new urban parameterizations in ORCHIDEE, we compare simulated grid cell averaged fluxes with observed fluxes measured at urban flux tower sites from the Urban-PLUMBER intercomparison (Lipson et al., 2022b). Following the Urban-PLUMBER protocol, the AU-Preston site in Melbourne, Australia, is used as a reference station for detailed evaluation (Lipson et al., 2024). In addition to evaluating the diurnal cycle at AU-Preston, we extend the analysis to 20 urban sites spanning a range of climates and urban fractions. For these sites, we analyse the distribution of daily maximum and minimum fluxes across summer and winter seasons to assess the robustness of model behaviour under contrasting energy balance conditions.

Model performance is quantified using two complementary metrics: bias and mean absolute error (MAE). The bias is defined as the mean difference between simulated and observed fluxes (sim − obs), providing information on systematic over- or underestimation. The MAE is defined as the mean of the absolute differences between simulated and observed fluxes, thereby measuring the overall magnitude of deviations regardless of sign. Bias is computed for sensible heat flux (H), latent heat flux (LE), and net radiation (R_net), and MAE is computed for sensible and latent heat fluxes. Ground heat flux (G) is also analysed in the model outputs, although no observational data are available for evaluation.

No additional statistical tests are applied, as the focus is on systematic model–observation differences across multiple years of flux tower data. By comparing results across climates and varying degrees of imperviousness, we assess the consistency of model improvements and identify the conditions under which the new scheme performs best.

4 Results

4.1 Energy budget

Even without the three-dimensional effects of buildings and the source of anthropogenic heat fluxes, urban areas significantly alter the surface energy budget compared to natural landscapes, typically modifying net radiation, increasing nocturnal sensible heat due to building material, and modifying the partitioning between sensible, latent and ground heat fluxes. To assess whether our model developments accurately capture these differences, we first compare simulated fluxes with observations from the AU-Preston station in Melbourne, Australia (Fig. 3), the reference site used in both Urban-PLUMBER intercomparisons (Best et al., 2015; Lipson et al., 2024).

Figure 3 Diurnal cycle of net radiation (R_net), sensible heat (H), latent heat (LE), and ground heat flux (G) during summer months (DJF – first row), and winter months (JJA – second row) of observational period August 2003–November 2004 over the station AU-Preston in Melbourne, Australia, the reference station for the Urban-PLUMBER intercomparison projects. The black line represents the observations at the station, while the coloured lines represent the simulations outputs. G was not measured directly; the observational estimate shown here was diagnosed as the surface energy-balance residual, $G=H+\mathit{LE}-R_{\mathrm{net}}$.

The simulations show small differences in net radiation in both summer and winter (Fig. 3). However, looking at upwelling shortwave (SW_up) and longwave (LW_up) components reveals compensating differences among the configurations. The grid-cell albedo increases from 0.110 in Baresoil and 0.131 in Lowveget to 0.141 in Urban0, Urban1, and Urban2. The albedo values discussed here correspond to the effective grid-cell albedo used to calculate net radiation, rather than to the albedo assigned to an individual PFT. The grid-cell albedo is calculated as the land-cover-fraction-weighted mean of the PFT-level albedos. Because of this increase, the urban configurations reflect more shortwave radiation than the benchmark simulations: Baresoil produces the lowest SWup, Lowveget intermediate values, and the three urban configurations the highest values. The urban configurations also modify the diurnal cycle of upwelling longwave radiation. Relative to Baresoil and Lowveget, they emit less longwave radiation around midday and more during the night, particularly in summer. This pattern indicates a reduced diurnal amplitude of surface radiative temperature and is consistent with the greater effective thermal inertia resulting from the changes in heat capacity and thermal conductivity. For the ground heat flux, urban simulations exhibit lower negative daytime values (more heat entering the column) and higher nighttime values (more heat released). More specifically, the Baresoil and Lowveget simulations produce overlapping curves, and the three urban configurations (Urban0, Urban1, and Urban2) likewise lie nearly on top of each other. These changes in LE, H and G, can be interpreted as a result of increasing heat storage capacity and thermal conductivity of urban surfaces. A higher C_p tends to reduce G, as it increases the amount of energy required to raise the temperature of the soil, leading to slower energy transfer into the ground. However, higher λ increases the rate at which heat is conducted into the ground, thereby increasing G. This increase in G consequently reduces the portion of energy available for sensible and latent heat fluxes at the surface, leading to lower energy partitioning towards atmospheric processes. Additionally, the three urban scenarios exhibit negligible differences, indicating that the representation of imperviousness has limited influence on the energy budget at AU-Preston.

When evaluated against the observations, all simulations slightly overestimate R_net around midday in summer. The radiative decomposition indicates that the urban albedo parameterization corrects most of the shortwave contribution to this bias; in the urban configurations, the remaining positive R_net bias is therefore mainly associated with an underestimation of upwelling longwave radiation around midday. However, the three urban scenarios (Urban0, Urban1, and Urban2) better reproduce the observed summer sensible heat flux, both at night (close to zero) and during the day. Despite this improvement, they still overestimate daytime sensible heat flux, which is likely related to the overestimation of R_net. These three urban scenarios still underestimate daytime sensible heat flux in winter. The comparable negative bias in H and positive bias in LE suggests that the main issue lies in the simulated partitioning of available energy between sensible and latent heat fluxes. Latent heat fluxes are reduced during the day but enhanced at night in urban simulations, with stronger imperviousness further reducing daytime values due to limited water availability. This diurnal contrast arises because, during the day, available energy is high but evaporation is limited by water availability, whereas at night, the release of heat stored in urban materials maintains relatively warm surface temperatures, allowing residual moisture to evaporate, which enhances the latent heat flux compared to natural surfaces. However, this effect is minor compared to other flux modifications. In winter, all simulations overestimate latent heat, though these discrepancies remain small relative to other energy fluxes. In summer, Baresoil and Lowveget underestimate the magnitude of daytime heat storage, as their simulated G values are less negative than the observational energy-balance residual. The urban configurations reproduce the observed daytime magnitude more closely, although their slightly more negative values indicate an overestimation of heat storage. At night, the urban configurations also agree more closely with the observational residual, but their more positive values suggest that nocturnal heat release is overestimated. This comparison remains qualitative because the colored curves represent the ground heat flux explicitly calculated by ORCHIDEE (Eq. 7), whereas the black curve is an energy-balance residual rather than a direct observation of G. The observational residual may additionally include contributions from unaccounted anthropogenic heat and net horizontal advective transport, as well as errors associated with measurement uncertainty.

To extend the analysis beyond a single location and assess whether the findings at AU-Preston are representative of broader conditions, we now examine 20 stations spanning different climates and urban covers. This allows us to evaluate whether the sign and magnitude of the model biases identified at AU-Preston (Fig. 3) are consistent across sites. Figure 4 presents the distribution of biases between maximum and minimum daily values of simulated and observed net radiation, sensible heat, and latent heat for winter and summer, two seasons with contrasting surface energy dynamics over the 20 urban sites. For each site, simulation, season, and flux, the daily simulated-minus-observed differences are calculated using only matching model–observation timestamps and are then averaged over all valid days. Each site therefore contributes one value to each boxplot, so that all sites are weighted equally regardless of the length of their observational record. Because the characteristic magnitudes differ among fluxes and between daytime maxima and nighttime minima, the absolute bias magnitudes should not be interpreted as a direct ranking of relative model performance across panels. The comparisons below focus primarily on differences among model configurations within each diagnostic. This approach highlights model performance during both daytime and nighttime phases of the diurnal cycle.

Figure 4 Distribution of mean bias in simulated daily maximum and minimum surface energy fluxes (sensible heat H, latent heat LE, and maximum net radiation R_net) during summer (first row) and winter (second row) months (DJF/JJA based on hemisphere) across the 20 urban sites. For each boxplot, the central horizontal line indicates the median, the lower and upper edges of the box indicate the 25th and 75th percentiles, respectively, and the whiskers extend to the most extreme values within 1.5 times the interquartile range (IQR). The dashed horizontal line indicates zero bias.

For maximum daily net radiation, the Urban simulations exhibit less bias across the stations, indicating an improvement brought by specifying the observed urban albedo (see Supplement for grid cell albedo differences across the stations and simulations). In contrast, the Baresoil and Lowveget simulations show biases with more variability, with an overestimation by Baresoil simulations in winter and an underestimation by Lowveget simulations in summer.

The different parameterizations used in the model configurations lead to the largest differences for the maximum daily sensible heat flux. During the summer, the Baresoil simulations have a more pronounced overestimation compared to the urban simulations, showing a median bias close to 77 W m⁻². In contrast, the Urban0 and Urban2 simulations have a bias distribution centered closer zero (respectively having a median bias of 14 and 18 W m⁻²), and 50 % of biases between 66 and −36 W m⁻². The substantial change in simulated sensible heat is due to increased thermal conductivity in urban simulations, as already seen in Fig. 3. However, many simulations do not match the maximum daily observed sensible heat as closely as for AU-Preston, although the biases are significantly reduced. The Urban2 simulation introduces only minor changes compared to the Urban0 simulation for sensible heat. Urban1, on the other hand, increases the bias compared to Urban0 and Urban2 during summer, but still brings improvement compared to Baresoil. Still, all simulations tend to overestimate sensible heat, the bias exceeding 150 W m⁻² at some stations. The largest summer sensible heat flux biases are common across the Baresoil, Lowveget, and Urban simulations. Additional sensitivity experiments indicate that they are partly related to underestimated ground heat uptake associated with the representation of soil thermal properties: prescribing higher thermal conductivity and heat capacity substantially reduces the biases at most affected sites, all of which have urban fractions below 50 %, highlighting the soil thermal parameterization as an important target for future ORCHIDEE development (Verhoef et al., 2026). In winter, maximum daily sensible heat is generally underestimated across the simulations. Baresoil has the least negative bias, although its bias distribution remains relatively broad, while Lowveget and the Urban simulations show similarly negative biases. This similarity suggests that the winter underestimation of sensible heat cannot be primarily attributed to the increased thermal conductivity introduced in the urban configurations. At the same time, all simulations show a positive bias in maximum latent heat flux, indicating that excessive partitioning of the available energy toward latent heat may contribute to the negative sensible heat bias. The absence of anthropogenic heat fluxes from the simulations may additionally contribute to the underestimation at some highly urbanized sites. For example, Urban-PLUMBER reports mean anthropogenic heat fluxes of 78.5 W m⁻² at UK-KingsCollege and 92.7 W m⁻² at KR-Jungnang.

Minimum sensible heat flux, which corresponds to the nocturnal sensible heat flux, is generally higher in cities than in natural areas due to the release of stored heat accumulated during the day and the persistence of heat transfer by convection at night. As expected, both in winter and summer, Baresoil and Lowveget simulations exhibit a negative bias, meaning they simulate nocturnal sensible heat that is too low. The Urban0, Urban1, and Urban2 simulations lead to improvements in the simulation of nocturnal sensible heat flux. Furthermore, the Urban1 simulation, which includes stronger imperviousness, further improves the representation of nocturnal sensible heat during both seasons. Still, Urban0, Urban1, and Urban2 show negative bias in minimum sensible heat fluxes in winter. This underestimation is unlikely to be primarily explained by latent heat flux, for which nighttime biases remain comparatively small. It may instead partly reflect an insufficient release of stored heat at night, as well as missing anthropogenic heat fluxes at highly urbanized sites.

Looking at maximum latent heat flux, Fig. 4 shows a clear seasonal contrast. In winter, biases are positive for all simulations, indicating an overestimation of daytime latent heat flux. In summer, biases are instead generally negative, indicating an underestimation of maximum latent heat flux, although the distributions remain broad and span both positive and negative values across stations. The summer underestimation may partly result from an insufficient representation of active vegetation and from the absence of irrigation in the model. This interpretation is supported by the Lowveget simulation, whose median summer bias is close to zero and which also shows a maximum sensible heat flux bias centered around zero. However, the distribution of Lowveget latent heat biases remains very broad across stations, suggesting that highly local factors such as vegetation characteristics, vegetation management, and irrigation, which are not explicitly represented by the model configuration, strongly influence latent heat flux. The Urban0 and Urban2 simulations show only limited differences compared with Baresoil. In contrast, Urban1, with stronger imperviousness and therefore lower water availability for evapotranspiration, further increases the summer underestimation while reducing the winter overestimation. Biases in minimum latent heat flux are comparatively small but consistently positive in both seasons across all simulations, indicating an overestimation of nocturnal latent heat flux.

Overall, the urban configurations improve the representation of the surface energy balance through both a better simulation of net radiation and changes in energy partitioning. In summer, the large overestimation of maximum sensible heat flux in Baresoil appears to arise from combined biases: net radiation is imperfectly represented, while both daytime heat storage and latent heat flux are underestimated, leaving too much of the available energy to be transferred as sensible heat. The Urban configurations improve net radiation and substantially improve the representation of daytime heat storage, which strongly reduces the sensible heat bias. Maximum latent heat flux, however, remains generally underestimated, suggesting that evapotranspiration, including the representation of local vegetation characteristics and irrigation, is an important remaining source of error. In contrast, the Urban configurations do not improve sensible heat in winter, they exhibit a systematic underestimation of maximum sensible heat, consistent with an overestimation of maximum latent heat flux. Despite this limitation, the urban configurations systematically improve nocturnal sensible heat in both seasons relative to Baresoil and Lowveget.

4.2 Mean Absolute Errors of latent and sensible heat fluxes over the 20 urban flux towers

Following the Urban-PLUMBER intercomparisons (Best et al., 2015; Lipson et al., 2024), we evaluate model performance using the mean absolute error (MAE), which quantifies the average deviation between simulated and observed fluxes. We compare the MAE of the Baresoil simulation (representing the original, non-urbanized version of ORCHIDEE) with that of the Urban1 simulation, which incorporates the new urban parameterizations based on findings from the previous sections. Figures 5 and 6 present the MAE for latent and sensible heat fluxes respectively across the 20 urban flux tower sites for both simulations. On average across all sites, the Urban1 simulation outperforms the Baresoil simulation for both latent heat (by 2.9 W m⁻²) and sensible heat (by 5.1 W m⁻²). At the individual station level, Urban1 improves the MAE for latent heat at all sites. For sensible heat, improvements are also seen at nearly all sites, except for CA-Sunset, US-WestPhoenix, UK-Swindon, and FI-Kumpula, where performance is slightly degraded. Detailed analysis suggests that these degradations arise from site-specific changes in energy partitioning. At FI-Kumpula and US-WestPhoenix, improvements in net radiation are not accompanied by sufficient partitioning toward latent heat and/or heat storage, resulting in increased sensible heat flux. Their urban fractions (46 % and 48 %, respectively) lie just below the 50 % threshold used to activate the prescribed urban thermal conductivity and heat capacity. Consequently, the urban configurations modify the radiative and aerodynamic properties at these sites without activating the corresponding increase in substrate thermal inertia, which may contribute to insufficient daytime heat storage and the degradation in sensible-heat-flux performance. The new urban developments significantly improve performance at previously poorly performing sites, such as GR-Heckor (−25 W m⁻²), Au-Preston (−14 W m⁻²), and FR-Capitole (−13 W m⁻²) for sensible heat, and SG-TelokKurau (−15 W m⁻²) for latent heat. However, substantial biases remain at some locations, such as SG-TelokKurau and UK-KingsCollege, where sensible heat remains underestimated in both the Baresoil and Urban1 simulations. Including anthropogenic heat fluxes in future model versions may help to address these biases. Figures 5 and 6 also display the mean MAE values for sensible and latent heat from the latest Urban-PLUMBER intercomparison (Lipson et al., 2024) at AU-Preston, representing ensemble model performance using site-specific information. The results show that ORCHIDEE, with the new urban developments, achieves performance comparable to Urban-PLUMBER models for sensible heat flux, whereas further improvements are needed to reach comparable performance for latent heat flux.

Figure 5 Mean Absolute Error (MAE) of latent heat fluxes (LE) over the 20 urban tower fluxes for the Baresoil and Urban1 simulations. The Urban-PLUMBER value (one black triangle) corresponds to the median Mean Absolute Error of latent heat flux from the “all – detailed” experiment shown in Fig. 6 of Lipson et al. (2024). The colours of the station's labels along the x-axis indicate the climate of station: green for Temperate, purple for Continental, orange for Dry, and blue for Tropical. The percentage in these labels indicates the urban fraction around the station. The dashed horizontal lines indicate the mean MAE across all stations, shown in orange for Baresoil and in purple for Urban1.

Figure 6 Same as Fig. 5, but for the sensible heat fluxes (H).

4.3 Annual water budget

In addition to altering the surface energy budget, urban areas also modify the surface water budget. Most notably, urbanization tends to increase surface runoff and decrease evapotranspiration, as reflected in the reduced latent heat fluxes discussed in the previous section. While 1D simulations are valuable for comparing simulated fluxes with tower-based flux observations from the Urban-PLUMBER dataset, there are no corresponding observations of soil moisture or runoff available at these stations. Therefore, we limit our analysis to a comparison of the mean annual water budget between the natural surface simulations (Lowveget and Baresoil) and the Urban simulations that include impervious surfaces (Urban1 and Urban2). Figure 7 presents the mean annual water budget across the 20 stations for the four simulations. In addition to climate, initial soil textures and the total amount of annual precipitation (visible in Fig. 7 as the sum of evapotranspiration, surface runoff, and subsurface runoff) and the precipitation characteristics (not visible in Fig. 7, e.g., short intense events vs. long-duration low-intensity rainfall, convective vs. stratiform) also influence differences in the annual water budget between stations. Because runoff and subsurface-runoff observations are unavailable at these sites, this analysis should be interpreted as a sensitivity assessment rather than a validation of hydrological performance. Evaluating the simulated runoff would require spatially distributed simulations over instrumented urban catchments and comparison with observed streamflow, which is beyond the scope of this study.

Figure 7 Mean annual partitioning of precipitation into evapotranspiration, surface runoff, and subsurface runoff for each site and simulation. Each stacked bar represents the three components as fractions of mean annual precipitation. Mean annual precipitation is indicated above each panel. Simulated years depend on station (12–16 years). Titles' colours indicate the climate of station: green for Temperate, purple for Continental, orange for Dry and blue for Tropical. Percentage in title indicates the urban coverage around the station.

Looking at the natural land cover simulations (Baresoil and Lowveget), most stations show a water budget dominated by evapotranspiration, except in Canada, where subsurface runoff dominates. This is likely due to energy limitations rather than water availability, which constrains evapotranspiration at that site. In the Urban2 simulation, evapotranspiration is generally not changed compared to the Baresoil simulation (+1.1 %), but slight increases are observed at some stations (e.g., JP-Yoyogi and FI-Torni), while a decrease is found at the tropical site SG-TelokKurau. The small increases may be linked to enhanced nocturnal latent heat flux, as discussed in the previous section, where redistributed energy from deeper soil column during the night can contribute to evaporation.

The Urban1 simulation, which includes a higher degree of imperviousness, consistently shows a reduction in evapotranspiration (−7.6 %) compared to the Baresoil simulation across all stations. However, the dominant changes in the annual water balance introduced by the Urban simulations are not related to evapotranspiration, but rather to surface and subsurface runoff. In the natural land cover simulations, subsurface runoff exceeds surface runoff at more than half of the stations (11 out of 20), while the remaining stations (9) show higher surface than subsurface runoff. The Urban1 simulation, which assumes higher imperviousness, reverses this ratio by significantly reducing subsurface runoff and increasing surface runoff. Only five stations under this scenario still show subsurface runoff exceeding surface runoff. In addition to altering the ratio, Urban1 doubles the annual surface runoff at most stations. The Urban2 simulation, which represents a more moderate reduction in saturated hydraulic conductivity, generally increases surface runoff and decreases subsurface runoff, although the magnitude of these changes is much smaller than in Urban1. However, the response is not consistent across all sites: Urban2 slightly increases subsurface runoff at KR-Ochang and SG-TelokKurau and decreases surface runoff at FI-Torni and FI-Kumpula. These site-dependent changes indicate that the annual water-budget response to a moderate conductivity reduction is not necessarily monotonic. Although saturated hydraulic conductivity is scaled linearly with imperviousness, the resulting water-budget response may be nonlinear because runoff generation and vertical water transfer depend on infiltration thresholds, rainfall intensity, antecedent soil moisture, soil texture, evapotranspiration, and changes in soil-water storage. These results suggest that representing saturated hydraulic conductivity as a function of imperviousness yields physically consistent outputs only when conductivity is strongly reduced, as in the Urban1 hypothesis, but not when the reduction is moderate, as in Urban2.

5 Discussion

This work represents a first step in integrating urban areas into the ORCHIDEE Land Surface Model. By introducing several key developments, this research has laid the foundation for more accurate urban representation within the model, thereby advancing our ability to simulate urban-climate interactions. In terms of surface energy budgets, the new developments produced increased night-time sensible heat flux, decreased day-time sensible heat flux, reduced latent heat flux, and increased ground heat flux with heat released at night, leading to higher night surface temperatures. This is consistent with the established urban surface-energy-balance framework (Loridan and Grimmond, 2012; Oke, 1988). These changes, together with the reduced biases relative to observations at urban flux-tower sites, confirm that the new urban scheme provides a more realistic representation of turbulent fluxes. The Urban-PLUMBER intercomparison similarly found that relatively simple slab schemes can perform well when vegetation, water availability, and imperviousness are adequately represented, indicating that model skill is not determined by geometric complexity alone (Lipson et al., 2024). In the water budget, the influence on total evapotranspiration and total runoff remains modest; however, the scheme substantially modifies the partitioning between surface and subsurface runoff. This shift is particularly important in urban environments, where enhanced surface runoff and reduced subsurface runoff are key drivers of rapid hydrological responses.

Despite the implemented urban developments, the simulations exhibit persistent seasonal biases. Most notably, during summer, urban simulations show an overestimation of daily sensible heat flux and an underestimated latent heat flux, indicating that too much of the available energy is partitioned into sensible rather than latent heat. This high-sensible-heat, low-latent-heat bias is a recognized limitation of urban land-surface models and has been associated with the omission or insufficient representation of urban vegetation and water availability (Best and Grimmond, 2015). Indeed, latent heat flux remains underestimated by many current urban models (Lipson et al., 2024). In winter, the daytime bias reverses, with an underestimated sensible heat and an overestimated latent heat. The absence of anthropogenic heat in the current scheme may contribute to the winter underestimation of sensible heat, as its representation has been identified as important for simulating wintertime urban energy budgets (Hertwig et al., 2020; Jin et al., 2021; Karsisto et al., 2016). Future work will test whether prescribing anthropogenic heat fluxes from observation-based gridded datasets can reduce this bias in kilometer-scale simulations. On the contrary, the winter overestimated latent heat is not consistently reported across other models (Lipson et al., 2024) and may therefore reflect ORCHIDEE-specific seasonal controls on evaporation and water availability. At night, sensible heat also remains underestimated in winter, while latent heat exhibits smaller positive biases in both seasons. The bias pattern displayed by this new urban representation therefore suggests that evaporation is insufficient in summer, when urban vegetation and irrigation are most active, but insufficiently constrained in winter and at night. A potential development for ORCHIDEE would therefore be to embed an active vegetated fraction, including seasonal phenology, directly within the urban PFT. This would allow urban vegetation, and potentially irrigation, to enhance transpiration during the summer while limiting vegetation-related evaporation during winter.

These biases should also be considered alongside the inherent uncertainties in flux-tower measurements, which often include some degree of energy-balance non-closure (Mauder et al., 2020). This context helps frame, though does not fully explain, the model–observation differences. Several avenues for improving the one-tile urban scheme developed here emerge thus as clear perspectives. First, the evaluation of ground heat flux (G) revealed strong sensitivity to parameter choices, specifically the thermal conductivity, highlighting G as a key, yet often underrepresented component of the urban surface energy balance. For instance, G was absent from the Urban-PLUMBER harmonized gap-filled datasets for 20 urban flux tower sites used in this study (Lipson et al., 2022b). Our current thermal conductivity value (3.24 W m⁻¹ K⁻¹), taken from the Noah-MP land-surface model (He et al., 2023), is substantially higher than values commonly used in CLM (Lawrence et al., 2019) for example, 0.767 W m⁻¹ K⁻¹ in its original urban implementation (Loridan and Grimmond, 2012) or 1.55 W m⁻¹ K⁻¹ in the later SURY configuration (Wouters et al., 2016), and may therefore require further refinement. Future work could focus on calibrating thermal conductivity and heat capacity based on observed or region-specific data. While no global database of urban thermal properties was identified during this study, using measured values where available would likely improve realism. Another planned improvement is to update thermal conductivity and heat capacity not only in grid cells where urban cover exceeds a threshold (here 50 %), but by applying a weighted average of natural and urban values based on the actual urban fraction.

A more refined representation of urban hydrology could also help address the above biases. Urban land surface models use a range of runoff strategies, from routing a fixed fraction of rainfall directly to runoff to more process-based formulations accounting for infiltration capacity, soil saturation, surface-water storage, and transfers between pervious and impervious tiles. The simpler formulations may be insensitive to rainfall intensity, antecedent soil moisture, and site characteristics (Jongen et al., 2024). The present study introduces a novel physically based representation of imperviousness by linking it to saturated hydraulic conductivity, thereby affecting infiltration, subsurface water transfer, and runoff generation. This development directly addresses the weak sensitivity of runoff to imperviousness found in 7 of the 18 models evaluated by Jongen et al. (2026). Nevertheless, Jongen et al. (2026) also found that 10 models omitted at least one major runoff-generation mechanism. Future developments should therefore represent surface ponding and interception storage, irrigation, drainage connectivity, and water transfers between urban and vegetated tiles. These processes should ultimately be evaluated against hydrological observations, because the present sensitivity analysis alone cannot establish whether the simulated water partitioning is realistic. Importantly, advancing the hydrological realism of the model should not be guided by latent heat performance alone, but by a comprehensive evaluation of the full water balance (Jongen et al., 2024). However, such evaluation is currently limited by a lack of observational data for key hydrological variables (e.g., runoff and soil moisture) at the urban flux tower sites used in this study (Jongen et al., 2026, 2024). This underscores the need for coordinated efforts to collect hydrological data in urban environments to support model development and validation.

The developments presented here were carried out to enable future simulations for which ORCHIDEE, now including an urban representation, is coupled with an atmospheric model to investigate climate–urban interactions. One key perspective of this work is to conduct high-resolution, convection-permitting simulations over Paris to investigate urban effects on temperature, precipitation, and catchment water balance. As ORCHIDEE is used at various spatial resolutions, from coarse to convection-permitting high resolution (< 3 km), the development undertaken here is a one-tile scheme which can fit both uses. As highlighted by Zhao et al. (2021), despite the impact of urban areas on climate, most state-of-the-art Earth system models and general circulation models participating in the CMIP do not account for urban regions. This is due to the limited impact of urban areas at coarse resolutions (100–200 km). Although constructing multi-model ensembles is considered best practice for providing robust and well-characterized climate projections, relying on non-urban projections for assessing climate-related risks in built environments is inadequate due to their bias regarding urban impacts. To bridge this gap, recent efforts have focused on two approaches: physically downscaling climate projections from one or two GCMs using regional climate models with urban land schemes or directly incorporating urban representation into GCMs. Our development aligns with both approaches. At convection-permitting high resolution, the representation of urban areas becomes even more critical, as cities are explicitly resolved at this scale, and an inaccurate urban energy balance can lead to incorrect atmospheric feedbacks. The CORDEX Flagship Pilot Study URB-RCC is currently investigating how different urban schemes influence simulations of urban temperature and precipitation. Results from the Urban-PLUMBER intercomparison (Lipson et al., 2024) have already shown that performance in resolving the surface energy budget is comparable between one-tile and canyon-type urban schemes. For this reason, the development of a one-tile scheme offers a robust and computationally efficient option for representing urban surfaces in convection-permitting simulations, while remaining scalable for large-scale or global applications. Still, questions related to operational actions such as urban greening, depaving, or planning strategies require insights from more complex schemes that resolve the street canyon, which our urban scheme in ORCHIDEE won't be able to provide.

6 Conclusion

This study presents a first step in implementing urban areas in the ORCHIDEE land surface model. By moving from a baresoil representation of cities to a one-tile urban scheme, we introduce key urban surface properties and propose a novel imperviousness parameterization that links saturated hydraulic conductivity to urban land cover. The new Urban1 scheme improves the simulation of surface energy fluxes across the 20 urban sites of the Urban-PLUMBER intercomparison project, as reflected by reduced MAE in both latent and sensible heat compared to the Baresoil version. Across all stations, mean MAE (mean absolute error) decreases from 30.3 to 27.3 W m⁻² for latent heat and from 41.9 to 36.9 W m⁻² for sensible heat. At sites with more than 50 % impervious cover, the improvements are even larger, with latent heat decreasing from 31.5 to 27.7 W m⁻² and sensible heat from 46.3 to 38.2 W m⁻². Performance at AU-Preston is closer to the benchmark range of the Urban-PLUMBER ensemble. For sensible heat, where Urban-PLUMBER mean benchmarks typically achieve MAE between 18.6 and 26.1 W m⁻², the original Baresoil representation had an MAE of 47.1 W m⁻², while the new urban scheme (Urban1) reduced it to 33.3 W m⁻². For latent heat, where the Urban-PLUMBER mean benchmarks are 18.5–32.9 W m⁻², the Baresoil scheme produced 34.5 W m⁻² compared to 33.7 W m⁻² with the new scheme. In addition, the updated representation modifies local hydrology as expected by increasing surface runoff and decreasing evapotranspiration. While these developments improve several aspects of urban energy and water balances, some limitations remain. In particular, the simulated reduction in subsurface runoff raises questions about the representation of urban recharge, which can be found to increase in observational studies (Bhaskar et al., 2016). The lack of urban runoff or soil moisture observations at the evaluated sites currently limits direct hydrological validation. Future work should address these gaps by running high-resolution 3D simulations at the urban catchment scale, enabling evaluation against streamflow observations. Additional developments, such as spatially and temporally varying urban parameters (e.g., from WUDAPT), interception of rainfall on urban structures, refinement of thermal properties, and the inclusion of anthropogenic heat fluxes, will further enhance the performances of this scheme. Together, these improvements will support the use of ORCHIDEE in studying urban climate processes and feedback from city to regional scales.

A perspective of this work is to conduct high-resolution ORCHIDEE simulations over large urban areas with a river routing scheme. These simulations would allow the evaluation of the urban module in which urban parameters such as albedo, building height, and imperviousness can vary both spatially and temporally. In the present 1D experiments, these parameters were prescribed for each site. However, the scheme developed here is capable of reading and interpolating urban data maps (e.g., LCZ maps or datasets describing imperviousness, albedo, and building heights), although this functionality was not demonstrated in the 1D framework. Using spatially distributed inputs would enable the model to take advantage of, and remain consistent with, comprehensive urban datasets such as the WUDAPT database (Ching et al., 2018). Including urban catchments within the simulation domain would also create opportunities to evaluate the model developments against streamflow measurements, which are generally more available than runoff or soil moisture observations. Further consideration could be given to whether additional processes, such as subsurface water lateral fluxes should be included in high-resolution climate simulations to assess urban water resources.