The experiment was carried out at the Jardin du Palais Universitaire in Strasbourg, France (48°35′04.4″N 7°45′48.2″E, 138 m altitude) from January 1^st, 2016 to December 31^th, 2016. The site consists in a vegetated place surrounded by some buildings. Its length and width are 220 and 90 m, respectively, and it is oriented at 112°. Meteorological variables and energy fluxes were measured by a 17 m height mast i.e, closed to the mean building height, thanks to several dedicated sensors. Samples were collected and stored every 15 min. A rotronic probe (HC2S3–Campbell) measured the air temperature and the relative humidity. A pyranometer CMP11 and a pyrgeometer CGR4 (Kipp and Zonen), directed downwards, measured the upwards solar and infrared radiations on the roof of a university building close to the park, respectively. A sonic anemometer (CSAT3 - Campbell) and a gas analyser (LI-7500A- LI-COR) measured the three wind speed components and the amount of water vapour in the air, respectively. LE and H were derived from Sonic data following the procedure proposed from Campbell. This place is a public garden with mainly grass, silver linden trees (Tilia tomentosa Moench), and some patches of bare soil (Figure 3).

Data were filtered in order to remove outliers by absolute deviation around the median method (Leys et al., 2013). This method uses median absolute deviation (MAD) following Eq. 1 (Huber, 2004). M A D = 1.4826 ∗ M i ( | x j − M j ( x j ) | ) (1) where M_j is the median of value at step t and the six observations around. M_i is the median of absolute difference between x_j and M_j. The value at step t is considered as an outlier and is removed from the dataset if it is outside interval spanning over the median plus/minus three MAD (Eq. 2). M j − 3 M A D < x i < M j + 3 M A D (2)

The area considered by the micrometeorological flux measurements, i.e., the so-called “footprint”, is changing with time and depends on multiple factors such as wind direction and speed, surface roughness, and thermal stability (Burba and Anderson 2010). In order to ensure that the areas considered by both measurements and TEB-Surfatm model were similar, the footprint analysis using the Neftel model (Neftel et al., 2008) was performed. Its aim was to determine the area contributing to at least 80% of the measured fluxes. The MAD filter on this subdataset was applied and called Footprint_MAD.

The annual trend of global radiation (Rg), in 2016, followed a typical pattern with minimum values during winter (mean daily values around 10 W/m²), increased during spring to reach a maximum in summer, and then decreased during autumn and winter. Mean daily Rg exhibited strong day-to-day variations, especially in summer (from around 100 W/m² to 350 W/m²) due to the interplay of cloudy and sunny periods (Figure 4A). At the daily scale, Rg increased from sunrise to reach its maximum (around 410 W/m² on average) at noon, and then decreased during the afternoon to reach its nocturnal value i.e., 0 W/m² (Figure 4B). The mean annual (±standard deviation) Rg was 130 ± 90 W/m², and 1^st and 3^rd quartiles were 40 W/m² and 200 W/m², respectively. Air temperature (Ta) and relative humidity (RH) were anti-correlated. At the daily scale, Ta increased from early morning to reach its maximum (around 18°C on average) in early afternoon and then decreased to its nocturnal value (around 12°C on average). RH followed the opposite trend by decreasing in the morning to around 60% on average in early afternoon before increasing to its maximum nocturnal value (around 90% on average) (Figure 4D). At the yearly scale, Ta followed the same pattern as Rg with minimum daily temperature around -4°C in winter and maximum around 28°C in summer. RH followed the converse annual trend and decreased from around 100% during winter to around 44% in summer before increasing during autumn (Figure 4C). As for Rg, RH and Ta exhibited large day-to-day variations i.e., drier and warmer days during sunny periods and wetter and colder days during cloudy periods. Mean annual Ta and RH were 12 ± 7°C (1^st and 3^rd quartiles equal to 6°C and 19°C) and 75 ± 13% (1^st and 3^rd quartiles equal to 65 and 85%), respectively. Over the whole year, rainfall accounted for 730 mm. These events occurred regularly, with 168 rainy days, with on average 9 ± 11 mm/day. Overall, it was dry in winter and wet in summer, with the wettest period from April to June. Indeed some particularly rainy days occurred during this period with daily rain accounting for up to 30 mm. The period from June to September was dry with only sparse rainy days (Figure 4G). The daily-mean wind speed varied between 0 m/s to 7 m/s and direction was mostly oriented to North-East (15–45°) and South (165–195°). The highest wind speeds occurred in the West-Southwest and East-Northeast directions (Figure 4H). Overall, the meteorological conditions were representative of the continental climate.

Daily means and mean half-hourly values of the measured net radiation (Rn) and sensible (H) and latent (LE) heat fluxes are represented in Figures 4E,F, respectively. The MAD filter (see Section 3.2) was applied on both H and LE data. At the daily scale, Rn, H, and LE followed the same dynamics: they remained stable during nighttime, increased during the morning to peak to their maximum at noon, and decreased to their nocturnal values during the afternoon. Maximum half hourly values were on average 300 W/m², 140 W/m² and 80 W/m² for Rn, LE, and H, respectively. During nighttime, Rn and H were both slightly negative while LE remained positive (Figure 4F). At the yearly scale, Rn followed the same seasonal and day-to-day dynamics as Rg. Mean daily Rn varied between around -25 W/m² and +25 W/m² during winter and between 75 and 190 W/m² during summer. Mean daily H did not exhibited seasonal evolution and ranged from -40 W/m² to 80 W/m², whereas mean daily LE increased from 0 W/m² during winter to 50–150 W/m² during summer and then decreased during autumn in link with the phenological development and physiological activity of the vegetation present in the site (Figure 4E).

In order to validate the TEB-Surfatm model, modelled energy fluxes were compared with the in-situ measurements which were filtered with the MAD method (see Section 3.2) i.e., by removing outliers.

Convective flux measurements were performed with the eddy-covariance method (see Section 3.1) i.e., a micrometeorological method providing integrated flux estimates on a certain area. This area, the so-called “footprint”, changes according to surface parameters (e.g., roughness or measurement height) and micrometeorological variables (e.g., wind direction and speed or atmospheric stability). Yet, TEB-Surfatm provides energy fluxes at three different scales: the vegetated area (including bare soil noted hereafter “Garden”), the road + garden area (hereafter “Canyon”), and the whole city including road, garden, walls, and roofs surfaces (hereafter “Town”). As a first approach, model outputs at these three different scales were compared to measurements in order to determine which scale was the most relevant for the comparison of modelled and measured fluxes. Four statistical indicators were used to assess these comparisons: the equation of the regression line, the Root Mean Square Error (RMSE), the bias, and the Mean Absolute Percentage Error (MAPE). The results are presented in Figure 5; Table 4.

Modelled and measured Rn exhibited a quite good agreement whatever the scale in model outputs considered i.e., garden (Figure 5A), canyon (Figure 5B), or town (Figure 5C) scale. However, the outputs at the canyon scale exhibited the best agreement with the measurements for both the regression line and the RMSE. Despite the best agreement with the measurements considering the bias and the MAPE statistical indicators for the town scale, this last also revealed the worst agreement considering the regression line and the RMSE. The Rn modelled at the garden scale exhibited the worst agreement with measurements considering bias and MAPE indicators (Table 4).

For the total convective heat fluxes (i.e., H + LE), although the fluxes estimated by the TEB-Surfatm model exhibited a quite good agreements with the measurements whatever the scale considered (Figures 5D–F), the best regression line was found for the fluxes modelled at the town scale which also exhibited the best RMSE but the worst bias and MAPE. On the contrary, total convective heat fluxes modelled at the garden scale exhibited the worst agreement considering regression line and RMSE. Those modelled at the canyon scale revealed the best agreement for both bias and MAPE (Table 4). Considering only latent heat fluxes, it was clear that the modelled LE at the town scale fitted the best with measurements since the best values for the four statistical indicators were found at this scale, and the four worst at the garden scale (Table 4). For the sensible heat flux, results were less clear: the town scale obtained the best regression but the worst values for the three other statistical indicators, the garden scale exhibited the worst regression but the best bias, and the canyon scale revealed the best values for RMSE and MAPE (Table 4).

Considering the whole results and statistical indicators, modelled fluxes at the canyon scale provided six best values for the statistical indicators, ten intermediates, and no worst values, while it was one best value, six intermediates, and nine worst values at the garden scale, and nine best, no intermediates, and seven worst values at the town scale (Table 4). These findings suggested that the canyon scale was overall the best modelling scale for the comparison with measurements. In addition, the fluxes were measured close to the mean building height. In this case it is a reasonable hypothesis to suppose that the fluxes mainly originated from the road and garden surfaces i.e., at the canyon scale. Therefore, fluxes are studied at the canyon scale in the following parts.

The partitioning of total convective heat fluxes in H and LE exhibited underestimation and overestimation of the measurements, respectively. Many issues could explain that finding such as the impact of the footprint from measurements or some missing or not well represented processes in the TEB-Surfatm model. These issues are explored in the following.

As stated previously, only the modelled fluxes at the canyon scale are considered. The aim of this part is to investigate the discrepancies between the modelled and measured fluxes. Hence, the residues, as the difference between measured and modelled values (i.e., positive when the observed values are larger than the predicted ones, and conversely), were estimated. The data were filtered with the MAD method (see Section 3.2; noted hereafter “MAD”), and also by including the filter based on the measurement footprint (noted hereafter “Footprint_MAD”). This latter consisted in keeping only the data for which at least 80% of the measured flux originated from maximum 200 m (i.e., the garden length) from the mast.

The residues were first analysed as a function of the deciles of measured energy fluxes. The Figures 6A,B show the boxplots of the residues according to the deciles of the measurements for H and LE, respectively. Typically, 1^st decile corresponds to low flux magnitude (i.e., during nighttime and/or winter) and 10^th decile corresponds to maximum flux magnitude (i.e., at noon and/or summer). Both results with MAD and Footprint_MAD filters are presented. Results of the pairwise comparisons are also shown: a letter shared between two or more boxplots indicated no significant differences, and two boxplots that do not share a common letter are significantly different.

The modelled and measured H fluxes presented a quite good agreement most of the time, residues boxplots being centred and closed to zero from 3^rd to 8^th measured H deciles (Figure 6A). Yet, considering the footprint for this measurement range did not improve the model accuracy since overall no significant differences where observed between residues with MAD filter and with Footprint_MAD filter for 3^rd to 7^th deciles of measured H. However, the model overestimated (i.e., negative values of the residues) H fluxes for low H values (i.e., 1^st and 2^nd deciles) (Figure 6A). In this case, applying the Footprint_MAD filter clearly and significantly improved the model/measures agreement. For the highest values of measured H (i.e., 10^th and to a less extend 9^th deciles), the model underestimated (i.e., positive values of the residues) the fluxes (Figure 6A). Applying the Footprint_MAD filter induced an even low but significant worst agreement with measurements (residues median, 1^st, and 3^rd quartiles equal to 46 W/m², 22 W/m², and 69 W/m², respectively) than with the MAD filter (residues median, 1^st, and 3^rd quartiles equal to 35 W/m², 12 W/m², and 56 W/m², respectively) (Figure 6A), suggesting that the measurement footprint was not the only reason for model overestimation. It will be explored in the following sections.

For the lowest LE values (i.e., 1^st decile of the measurements), the model strongly overestimated (i.e., negative values of the residues) the measurements. When LE fluxes were quite low (i.e., 2^nd and 3^rd deciles of measurements), the model agreed well with the measurements: residues were closed to zero whatever the filter (MAD or Footprint_MAD) considered. From 4^th to 10^th deciles of the measurements, the model globally overestimated the measured LE fluxes and its accuracy decreased when measured LE increased (i.e., the interquartile range (IQR) of the residues boxplots increased) (Figure 6B). Moreover, although applying the Footprint_MAD filter clearly improved the model/measures comparison for the highest values of LE (i.e., 9^th and 10^th deciles of measurements), it led to worst agreement for intermediates values of LE (i.e., 4^th, 5^th, and 6^th deciles of measurements) (Figure 6B).

Hence, the ability of the TEB-Surfatm model to estimate H and LE fluxes varied according to the flux magnitude. The model reproduced well H fluxes excepted for very high H fluxes and overestimated overall LE fluxes. The discrepancies partially originated from the measurement footprint: the area contributing to the measured flux could be larger than the area considered by the TEB-Surfatm model. However, the measurement footprint did not explain all the differences between measured and modelled fluxes. Other issues, such as biophysical processes not well taken into account could not be discarded. Correlations between residues and measured meteorological variables (global radiation, air humidity, wind direction and speed, soil temperatures, and soil water content) were tested (data not shown) but no significant correlation were established. In the following, the temporal (i.e., daily and seasonal) and spatial variabilities of the residues were explored to explain the differences between model and measured fluxes.

The two-hours residues boxplots for H and LE fluxes with MAD and Footprint_MAD filters are presented for spring (Figures 7A,B), summer (Figures 7C,D), and autumn (Figures 7E,F) seasons.

For H fluxes, the IQR of the residues exhibited a clear daily variation. It was weak during nighttime, increased during the morning, peaked to its maximum at noon, and decreased during the afternoon to the minimum IQR. This trend occurred whatever the season considered (Figures 7A,C,E). Hence, the absolute error between measured and modelled H fluxes varied throughout the day, but it must be kept in mind that the flux magnitude also varied with a maximum near noon (see Section 4.1 and Figure 4F). Hence, in terms of relative error, the performance of the TEB-Surfatm model was overall constant at the daily scale. Most often a slightly better agreement (median closer to zero, smaller IQR) is observed with the Footprint_MAD filter. Yet, even if it was quite weak, it must be noted that the TEB-Surfatm model slightly underestimated (i.e., positive residues) the H fluxes during spring (Figure 7A) and slightly overestimated (i.e., negative residues) the H fluxes during summer (Figure 7C) and autumn (Figure 7E).

As for the LE fluxes, the IQR of the residues exhibited a similar and clear daily variation for all seasons (Figures 7B,D,F) i.e., the absolute error of the model increased with the magnitude of the LE flux. Contrary to H fluxes, the model systematically overestimated the LE flux measurements at the two-hours time scale. Nevertheless, considering the measurement footprint improved the ability of the model to reproduce the measurements. With the Footprint_MAD filter, the residues boxplots were systematically closer to zero, but no dependence on the season occurred. Yet, it overall did not change the IQR of the residues boxplots.

The footprint depends notably on wind direction, therefore, the ‘measurement area’ is not centred on the mast and can largely differ from the ‘modelling area’. In order to investigate that issue, a method of bivariate interpolation and smooth surface fitting was used to present the residues values at points distributed in a x-y plane determined by distances and directions from the mast (Akima, 1978, 1996). These results are presented in Figure 8. Each pixel delimits the upwind distance from the mast for which 80% of the measured flux originated from. For each pixel, estimated residues was sorted in three categories for H (Figure 8A) and LE (Figure 8B).

This representation allows an interpretation of the differences between the model and measurements according to distance and direction from the mast i.e., the footprint of the measurements.

The H fluxes were underestimated (i.e., residues larger than 30 W/m²) by the model when the footprint was very closed to the mast, extended at 50 m south from the mast, and 80–90 m west-south-west from the mast. Model overestimation (i.e., residues lower than -30 W/m²) mainly occurred when the footprint extended to up to 100 m from the southern, western, and northern directions. However, when the measurements footprint was in the canyon considered by the TEB-Surfatm model, neither over- nor underestimation was observed (i.e., residues between -30 W/m² and 30 W/m²) (Figure 8A).

The modelled LE fluxes were underestimated (i.e., residues larger than 50 W/m²) mainly when the measurement footprint originated from south-south-west to north-west directions, and overestimated (i.e., residues lower than -50 W/m²) when it originated from north-east to east directions from the mast. The modelled LE fluxes exhibited a good agreement with measurements (i.e., residues between -50 W/m² and 50 W/m²) when the measurements footprint was in the canyon, except for the area close to the mast in the east direction (Figure 8B).

Hence, the variability of the model accuracy depended on the measurement footprint in link with the heterogeneity of the site characteristics. While the measurement footprint fitted with the area characteristics considered by the TEB-Surfatm model, the modelled H and LE fluxes presented a good agreement with the measurements. It proved the ability of the coupled model to estimate the energy fluxes in urban area, accounting for impervious and vegetated surfaces. The underestimation and overestimation in the western area from the mast for H and LE fluxes, respectively, could be attributed to a more important fraction of impervious surface than what it was considered in the model. Indeed, this zone consisted in mainly a built-up area (Figure 3). Therefore, the evapotranspiration is very weak and most of the energy received by the surface would be dissipated through H flux compared to a vegetated area (Gunawardena et al., 2017; Chrysoulakis et al., 2018). The model overestimation of LE fluxes when the footprint originated from north-north-east and east directions could also be attributed to a difference between the site and modelled characteristics. Indeed, the amount of trees in these areas is larger than the amount considered in the modelled area (Figure 3). Since the maximum stomatal conductance (g_max) is three times larger for grass than for tree (Table 3), the larger amount of grass in the model can lead to model overestimation of the LE fluxes in these areas. That issue was also observed in previous studies (e.g., Wang et al., 2012; Zou et al., 2021). For instance, Litvak et al. (2014) reported that the evapotranspiration of a grass plot was between +50% and +100% larger than a mixed grass and trees plot in summer.