Evaluation of the Soil, Vegetation, and Snow (SVS) Land Surface Model for the Simulation of Surface Energy Fluxes and Soil Moisture under Snow-Free Conditions

Abstract

The recently developed Soil, Vegetation, and Snow (SVS) land surface model is being progressively implemented at Environment and Climate Change Canada (ECCC) for operational numerical weather and hydrological predictions. The objective of this study is to evaluate the ability of SVS, in offline point-scale mode and under snow-free conditions, to simulate the surface heat fluxes and soil moisture when compared to flux tower observations and simulations from the Canadian Land Surface Scheme (CLASS), used here as a benchmark model. To do this, we performed point-scale simulations of between 4 and 12 years of data records at six selected sites of the FLUXNET network under arid, Mediterranean and tropical climates. At all sites, SVS shows realistic simulations of latent heat flux, sensible heat flux and net radiation. Soil heat flux is reasonably well simulated for the arid sites and one Mediterranean site and poorly simulated for the tropical sites. On the other hand, surface soil moisture was reasonably well simulated at the arid and Mediterranean sites and poorly simulated at the tropical sites. SVS performance was comparable to CLASS not only for energy fluxes and soil moisture, but also for more specific processes such as evapotranspiration and water balance.

1. Introduction

Land surface models (LSMs) simulate energy and water exchanges between the different components of a land system. LSMs are typically integrated into numerical weather prediction and climate projection schemes in order to provide the required surface boundary conditions [1,2]. They are also a component of land–surface data assimilation systems [3,4,5], as well as hydrological prediction systems [6].

LSMs for use in weather and climate modeling have been the subject of intense development over the last three decades [7,8]. They consolidate knowledge acquired from micrometeorology, soil sciences, as well as field and laboratory research in plant physiology. The complexity inherent to these models makes them remarkable in their capacity to capture energy and moisture exchanges within the soil, vegetation and snow components, as well as hydrological processes at the catchment scale [2]. Substantial progress associated with LSM development has also resulted from model inter-comparison studies [9,10,11] that highlighted the interest in running LSMs decoupled from their host atmospheric model, since offline simulations allow computationally inexpensive research and model improvement [12,13]. A good example is the Canadian Land Surface Scheme (CLASS) [14,15] that has been successfully assessed (e.g., [16,17,18,19]) and widely applied over the North American boreal zone (e.g., [20,21,22,23,24,25]) and other environments (e.g., [26,27,28,29,30]), and more recently in a hydrological context (e.g., [31,32,33]).

According to Pitman [7], the two key equations representing the land surface processes are the energy and water balances. The first one, applied to a surface, is defined as follows:

$$R_n=H+LE+G,$$

(1)

where $R_n$ is net radiation; H and $LE$ are the sensible and latent heat fluxes, respectively; and G is the soil heat flux. The water balance (under snow-free conditions) is defined as:

$$P=ET+R_{surf}+R_{base}+\mathsf{\Delta }S,$$

(2)

where P is precipitation, $ET$ is evapotranspiration, $R_{surf}$ is the overland flow (possibly including subsurface lateral flow), $R_{base}$ is the drainage (vertical flow from the bottom of the soil column) and $\mathsf{\Delta }S$ is the variation in the water stored in the surface scheme (i.e., in the vegetation canopy, snowpack or soil column). Note that $ET$ is linked to the surface energy balance by the latent heat flux term, $LE=\lambda ET$, where $\lambda$ is the latent heat of vaporization. LSM assessment is typically done in the context of the evaluation of all terms of Equation (1). [9,34,35,36] and/or Equation (2) [37,38,39].

Environment and Climate Change Canada (ECCC) currently uses a modified version of the Interaction between Soil, Biosphere and Atmosphere (ISBA) LSM to produce operational weather forecasts [40,41,42]. ISBA provides surface fluxes of momentum, water and energy to the Global Environmental Multiscale (GEM) atmospheric model. To address particular limitations with their operational version of ISBA [40,41], ECCC developed the Soil, Vegetation and Snow (SVS) LSM. As documented by Alavi et al. [43] and Husain et al. [44], SVS introduces the following improvements compared to ISBA: a new tiling approach, an improved parameterization of the vegetation thermal coefficient, the inclusion (as an option) of photosynthesis processes, a new formulation of land surface albedo and emissivity, a root density function that depends on the vegetation type, an update of the snow melt process and multilayer vertical water transport in the soil.

In a weather prediction context, SVS improved ISBA offline simulations of soil moisture and brightness temperature during the warm season over the North American domain [43,44]. In a hydrologic context, SVS has been shown to produce satisfactory runoff simulations within the large and partly ungauged Lake Ontario basin [6], and in hindcasts of the June 2013 flood in Alberta (Canada) [45]. More recently, Maheu et al. [46] evaluated the point scale soil water transport component of SVS by using precipitation and latent heat fluxes as input data. The model component performed well when simulating soil moisture in Mediterranean and temperate climates of the United States. While previous studies have evaluated the general performance of SVS and contrasted it to ISBA, no detailed evaluation of SVS outside of the North American domain has been performed to date. Moreover, none of the studies focused on evaluating the ability of SVS to simulate surface energy fluxes. Thus, there is a clear need for additional evaluations of SVS focusing on the model’s potential to simulate vertical exchanges of heat and water for various climates and biome types. More particularly, we are interested in the evaluation of soil and vegetation components, but not the snow component, which deserves particular attention and would be the subject of another study.

The objective of this study was thus to evaluate the ability of SVS to simulate surface energy fluxes and soil moisture in snow-free conditions for various climate and biome types. SVS was assessed in its offline point-scale configuration by comparing its outputs to energy fluxes and near-surface measurements. It was also compared to simulations from the widely implemented CLASS model, which was used as a benchmark model.

2. Methodology

2.1. Study Sites

FLUXNET is a network of micrometeorological measurements. It includes more than 800 active and historic flux measurement sites spread across most of the world’s climatic regions and representative biomes [47]. Six FLUXNET sites were selected from a total of 212 sites in the FLUXNET2015 dataset updated in November 2016. The sites were filtered progressively using the following four criteria: (1) yearlong snow-free conditions (50 sites); (2) data availability of surface energy fluxes and volumetric soil moisture content measurements over multiple years (12 sites); (3) contrasting climate types; and (4) biomes (6 sites, Figure 1).

Table 1. Main characteristics of the study sites.

Characteristic	AU-ASM	US-Wkg	US-Ton	IT-Noe	MY-PSO	GF-Guy
Latitude (${}^{\circ }$)	−22.28	31.74	38.43	40.61	2.97	5.28
Longitude (${}^{\circ }$)	133.20	−109.94	−120.97	8.15	102.31	−52.92
Altitude (m)	606	1530	177	25	150	48
Climate classification ${}^1$	Bsh	BSk	Csa	Csa	Af	Af
Biome type ${}^2$	ENF	GRA	WSA	CSH	EBF	EBF
MAP (mm) ${}^3$	306	407	559	588	1804	3041
MAT (${}^{\circ }$C)	20.0	15.6	15.8	15.9	25.3	25.7
Data period of simulation	2011–2014	2011–2014	2003–2014	2004–2014	2003–2009	2004–2014

${}^1$ Köppen–Geiger climate classification: hot semiarid (Bsh) and cold semiarid (BSk). ${}^2$ Biome types according to IGBP land cover classification: evergreen needleleaf forest (ENF), grassland (GRA), woody savanna (WSA), closed shrubland (CSH) and evergreen broadleaf forest (EBF). ${}^3$ Mean annual precipitation (MAP) and mean annual air temperature (MAT).

summarizes the main characteristics of the study sites, which were grouped, according to the Köppen–Geiger climate classification, as arid, Mediterranean and tropical climates. On the other hand, biome type was grouped according to the International Geosphere-Biosphere Programme (IGPB) as evergreen needleleaf forest (ENF), grassland (GRA), woody savanna (WSA), closed shrubland (CSH) and evergreen broadleaf forest (EBF).

The sites show mean annual precipitation, ranging from 306 mm at the arid AU-ASM site to 3041 mm at the tropical GF-Guy site (Table 1). The mean annual air temperatures range from 15.6 to 25.7 ${}^{\circ }$C (Table 1) The AU-ASM site is covered with a mulga savanna woodland at the southern end of the North Australian tropical transect [48]. At this site, 70% of its annual precipitation occurs during the summer months (December–February) and 85% falls throughout the monsoon season (November–April). The mean monthly temperature can range from 13 ${}^{\circ }$C in June and July to 33 ${}^{\circ }$C in January (Figure 2a). The other arid site, US-Wkg, exhibits a drastic seasonality in precipitation patterns (Figure 2b), with 70% of the annual precipitation taking place during the summer months from July to September because of the annual cycle of the North American monsoon [49]. The land is covered mainly by a diverse mosaic of native grassland species replaced progressively by Lehman lovegrass starting in 2006 [50]. The US-Ton Mediterranean site has a mean annual precipitation of 559 mm, with 80% of it concentrated from November to April and with dry and warm conditions during the summer months from July to September (Figure 2c). Land cover is classified as an oak savanna woodland (dormant in winter and active in spring and summer) interspersed with grassland and forbs (mainly active during the rainy winter periods) [51]. The other Mediterranean site is in Italy, IT-Noe, and displays a similar behavior compared to the US-Ton site in terms of annual temperatures. Precipitation mainly occurs during Spring and high water deficit is observed from May through September. Intense storms provide significant runoff and low water storage during fall [52]. The main species of vegetation are juniper and lentisk with variable-sized patches of bare ground [52]. The MY-PSO tropical site has a mean annual precipitation of 1804 mm, which is less than other regions of Peninsular Malaysia [53]. Rainfall peaks in March and April and from October to December (Figure 2e). Monthly air temperature does not show seasonality and remains uniform around 27 ${}^{\circ }$C (Figure 2e). The area is covered by a mixed dipterocard forest with a continuous canopy height of around 35 m and some emergent trees surpassing 45 m [54]. The other tropical site, GF-Guy, presents intense precipitation during December–February and April–July due to the north–south movement of the Inter-Tropical Convergence Zone [55]. This site is located in a tropical wet forest with a mean tree height of 35 m and emerging trees exceeding 40 m [55].

2.2. Data

Micrometeorological forcing and validation data were both retrieved from FLUXNET2015 dataset (www.fluxnet.fluxdata.org) at 30-min time steps. The forcing data consisted of seven atmospheric variables needed to pilot both SVS and CLASS: air temperature ( ${}^{\circ }$C), precipitation (kg m${}^{-2}$ s${}^{-1}$), downward shortwave radiation (W m${}^{-2}$), downward longwave radiation (W m${}^{-2}$), wind speed (m s${}^{-1}$), surface pressure (Pa) and specific humidity (kg kg${}^{-1}$). Specific humidity was estimated as a function of relative humidity and air temperature following Gill [56]. FLUXNET2015 data are gap-filled with a combination of autocorrelation gap-filling techniques [57] and downscaled ERA-interim data [58].

The surface energy fluxes to be evaluated are the latent heat flux ($LE$), sensible heat flux (H), net radiation ($R_n$), and soil heat flux (G). They were gap-filled using autocorrelation filling techniques [57]. The turbulent fluxes ($LE$ and H) were measured using the eddy-covariance technique [59]. They were also adjusted in the FLUXNET2015 dataset by an energy balance closure correcting factor ( ($R_n-G$)/($H+LE$) ). The factor avoid changes in the heat storage in the ecosystem. It was calculated at 30-min time intervals and applied only when energy closure fell outside 1.5 times the interquartile interval and when all the terms were available. Following similar studies by Blyth et al. [12] and Alves et al. [30], a daily energy balance closure filter was applied in order to guarantee a proper observed surface energy balance closure. The daily filter consisted of discarding days for which the error in energy balance closure exceeded 10%. The percentage of the data accepted is the following: AU-ASM (43.2%), US-Wkg (17.3%), US-Ton (11.2%), IT-Noe (29.4%), MY-PSO (91.6%), and GF-Guy was not filtered because observations of G were not available.

Observed volumetric soil water content corresponds to the topsoil layer, typically between 5 and 10 cm. It was also gap-filled using autocorrelation filling techniques [57]. Soil moisture observations were derived from electrical property of the soil. More exactly, theta ML2-X probes (Delta-T Devices) were deployed at US-Ton [60], Stevens Water Hydra Probe sensors at US-Wkg [61]; reflectometry probes at AU-ASM (CS605 and TDR100 Campbell Scientific [62]), at IT-Noe (CS616-L Campbell Scientific [52]), at MY-PSO (CS615 and CS616, Campbell Scientific [54]), and at GF-Guy (TRIME FM3; Imko, Ettlingen [55]).

2.3. Land Surface Models

2.3.1. SVS

SVS is a land surface model recently developed at ECCC [43,44]. In SVS, each grid cell may be divided into the following four tiling components: (1) bare ground; (2) vegetation (low and high); (3) snow over bare ground and low vegetation; and (4) snow under high vegetation. In the absence of snow, only Components (1) and (2) are taken into account. The model considers a single-layer approach for the vegetation canopy, a one-layer snowpack and a user-defined number of soil layers, with user-defined depth (by default set to seven layers) to compute soil moisture. The force-restore method is used to solve the energy budgets of each surface tile. Schematic diagrams of the hydrological regime and the tiling implementation in SVS soil part can be found in the works of Alavi et al. [43] and Husain et al. [44].

The energy budgets for the different tiles are obtained using the frce-restore approach. In the presence of different vegetation types, in a single grid cell, a single energy budget is calculated using mean vegetation parameters (such as albedo and emissivity) based on a weighted mean of the parameters of each individual type of vegetation present in the grid cell, and their respective fractional coverage. The net radiation for each type of vegetation is estimated taking into account total incoming solar and infrared radiation, as well as vegetation albedo and emissivity. The sensible and latent heat fluxes over vegetation are modeled using the bulk aerodynamic approach [43]. The soil heat flux was computed as the residual of the net radiation and the sensible and latent heat fluxes.

As described by Husain et al. [44], SVS distinguishes among 21 types of vegetation. In each grid cell, a varying combination of these vegetation types can be present. Vegetation parameters, such as leaf area index, roughness length, heat capacity of vegetation, albedo, stomatal resistance and root depth, are assigned separately for each vegetation-type and are specified using look-up tables. The canopy interception capacity is computed as a function of the leaf area index. If the canopy surface is covered with a film of intercepted rainfall (or if the vapor flux is toward the canopy), evaporation (condensation) takes place following the bulk aerodynamic approach. If the canopy is dry, the transpiration from leaves is computed by introducing the stomatal resistance to the evaporation. To compute the stomatal resistance, SVS has the option to use either a Jarvis-type approach [44] or a photosynthesis module that follows a biochemical approach [63,64,65]. In this study, the photosynthesis module was used.

The evolution of volumetric soil moisture in different soil layers is calculated using the diffusion equation. Liquid water flow is calculated using the Darcian equation for one-dimensional fluid flow along with the hydraulic conductivity and soil water suction following Clapp and Hornberger [66]. The saturated volumetric water content, wilting point volumetric water content and volumetric water content at field capacity are calculated based on the percentage of sand and clay contents in the soil according to Boone et al. [67] and Soulis et al. [68]. The saturated hydraulic conductivity and saturated soil water potential are calculated from the percentage of sand [69]. The lateral flow in the downhill direction is estimated through a parameterization of Richard’s equation, which estimates lateral flow rate at the seepage face of a hillslope. Lateral flow only occurs when soil water content exceeds water content at field capacity for a sloping element [68]. The water balance for the entire soil column is evaluated as a residual of precipitation, evaporation, water uptake by roots, runoff, lateral flow, baseflow and change in the soil water liquid content.

2.3.2. CLASS

CLASS was initially conceived in response to the perceived need for more soil thermal and moisture layers and for a separate treatment of vegetation canopy within the Canadian Global Circulation Model [14,15]. It has been under continuous development, as evidenced by Wu et al. [70], von Salzen et al. [71], Bartlett et al. [21], and Verseghy [72]. Each grid cell in CLASS may be divided into subareas of: (1) bare soil; (2) vegetation over the soil; (3) vegetation over snow; and (4) snow over soil. The vegetation canopy is treated as a single layer, same for the snowpack. The soil column has the flexibility to handle any number of soil layers to whatever depth the user specifies. The energy fluxes are computed based on the heat conservation equation described by Verseghy et al. [15]. The prognostic evolution of soil water content and the vertical water flow are given by the one-dimensional Darcy and Richards equations, respectively (see Verseghy [14] for a more detailed description.)

2.3.3. Differences between SVS and CLASS

As shown in

Table 2. Summary of SVS and CLASS main physical processes.

Process	SVS	CLASS
soil and vegetation turbulent fluxes	bulk aerodynamic approach [44]	bulk aerodynamic approach [14,15]
soil heat flux	residual of net radiation and turbulent fluxes	linear combination of soil temperatures [14]
surface (skin) and mean vegetation temperature	force-restore method [74,75]	heat diffusion equation [15]
transmissivity	Beer’s law [76]	Beer’s law [15]
canopy interception capacity	proportional to LAI [77]	proportional to LAI [15]
stomatal resistance	1. empirical parameterization [44,78,79] 2. Photosynthesis module (biochemical approach) [44,63,64,65]	empirical parameterization [15]
surface (skin) and mean soil temperature	force-restore method [74,75]	heat conservation equation [14]
soil moisture	moisture conservation equation [43,44]	moisture diffusion equation [14]
soil moisture flux	1-D Darcy’s law [43]	1-D Darcy’s law [14]
saturated hydraulic conductivity, saturated soil water suction	empirical equation based on soil texture [69]	empirical equation based on soil texture [69,80]
wilting point; field capacity; saturated values of volumetric water content, hydraulic conductivity, and soil water suction	parameterization from soil texture [67,69]	parameterization from soil texture [68,80]
relationship between soil water content, soil water potential, and relative permeability	empirical equation [43,66]	empirical equation [15,66]
vertical water flow	1-D Richards equation [81]	1-D Richards equation [14]
lateral water flow	1-D Richards equation [68,81]	1-D Richards equation [68,81]

, SVS and CLASS rely on similar parameterizations to represent several physical processes. However, there are two distinctive differences. One is associated with stomatal resistance, where the photosynthesis module is used in SVS, while an empirical parameterization is used in CLASS. Note that, although not in the version used in this study, CLASS can be coupled to a similar photosynthesis module [73] and can also have access to a climate vegetation model [71]. The other difference resides in the method used for computing the surface (skin) temperature and the mean temperature of the vegetation and soil components. A one-layer force-restore method is employed in SVS, while CLASS uses a multi-layer heat diffusion equation. Note that a multi-layer Richard’s equation framework for soil moisture is used in SVS and CLASS. The force-restore approach is still used in SVS for its computational efficiency and parsimony (less state variable) in the context of operational forecasting.

2.4. Experimental Setup

In this study, experiments were performed at point-scale and in offline mode with the MESH modeling platform (Modélisation Environnementale communautaire—Surface Hydrology), which is ECCC’s community environmental modeling system [31]. MESH allows different LSMs and hydrologic routing models to coexist within the same modeling framework. Although the final goal of MESH is to forecast water discharge at the watershed scale, it is also possible to utilize the stand-alone offline versions of SVS and CLASS, which are the two current LSMs integrated into MESH. A technical description of MESH can be found in Pietroniro et al. [31] and at www.wiki.usask.ca/display/MESH. MESH including the SVS and CLASS codes is available at www.wiki.usask.ca/display/MESH/Interim+releases.

Three time steps were specified in the models: (1) the input time step was set as 30-min for each model; (2) the model time step was defined as 5-min for SVS and 30-min for CLASS (a test carried out with CLASS for a time step of 5-min revealed almost the same results as those obtained at 30-min); and (3) the output time step was fixed as 30-min for both models.

The soil was divided into seven layers with depths from the top to the bottom equal to 10, 35, 85, 160, 235, 310 and 410 cm. Soil texture was assumed constant in the vertical direction since information at multiple depths was not available in the literature (

Table 3. Parameters for SVS and CLASS models at each FLUXNET site used in this study.

Parameters ${}^1$	AU-ASM		US-Wkg	US-Ton		IT-Noe	MY-PSO		GF-Guy
Vegetation type ${}^2$	ENF	LG	LG	EBF	LG	EBS	TEBF	LG	TEBF
Frac. coverage (% ) ${}^3$	30.0	45.0	50.0	40.0	30.0	75.0	70.0	15.0	80
Root depth (m)	1.0	1.2	1.2	5.0	1.2	0.2	5.0	1.2	5.0
Roughness length (m)	1.5	0.08	0.08	3.5	0.08	0.05	3.0	0.08	3.0
Measurement height (m)	11.7		6.0	23.0		2.0	52.0		52.0
Soil texture	sandy loam		sandy loam	silty loam		clay loam	loam		sandy clay
sand (%)	74		55	48		30	35		48
clay (%)	15		25	10		30	18		43
orgm (%) ${}^4$	1		0	0		0	10		10

${}^1$ References for the parameter values: AU-ASM: Cleverly et al. [62]; US-Wkg: Scott et al. [82]; US-Ton: Baldocchi et al. [51]; IT-Noe: Marras et al. [52]; MY-PSO: Kosugi et al. [54] and Yamashita et al. [83]; GF-Guy: Bonal et al. [55]. ${}^2$ Vegetation type for SVS and CLASS: ENF, evergreen needleleaf forest; LG, long-grass; EBF, evergreen broadleaf forest; EBS, evergreen broadleaf shrub; TEBF, tropical evergreen broadleaf forest. ${}^3$ Note: The remainder non-vegetative percentage corresponds to bare ground. ${}^4$ orgm, organic matter (only used in CLASS).

).

The soil texture parameters used at each site (Table 3) are provided to both models, except for the organic matter (orgm), which is used only by CLASS.

Vegetation associated parameters are defined in the look-up tables of each model. As shown in Table 3, both models account for the same vegetation type. For simplicity shown treating lateral water fluxes, the grid was considered perfectly flat at all sites.

Both models were set up with the same initial soil moisture and temperature conditions. Soil temperature for each layer was initialized with the yearly mean temperature. Soil moisture for each layer was initialized at 50% of saturation. To account for an equilibrium in the prognostic variables at each site (spin up), the first year of meteorological input was replicated ten times.

2.5. Modeling Performance

The performance of SVS and CLASS simulations was evaluated using the following metrics: root mean square error (RMSE), which measures the spread of the residuals between simulation and observation; percent bias (Pbias), which quantifies the average tendency of the simulated value to be larger or smaller than its observation counterpart; the Pearson correlation coefficient, R, which measures the degree of correlation between simulation and observation; and the Nash–Sutcliffe efficiency (NSE) for which the optimal value is 1 and a value inferior to 0 indicates that the average of the observation offers a better predictor than the model.

3. Result and Discussion

3.1. Energy Fluxes

For a precise evaluation of energy fluxes, the mean annual cycles of latent heat flux ($LE$), sensible heat flux (H), net radiation ($R_n$) and soil heat flux (G) are plotted in Figure 3. The metrics of the same fluxes were computed at a 30-min sampling interval and are summarized in

Table 4. Scores for the latent heat flux ($LE$), sensible heat flux (H), net radiation ($R_n$) and soil heat flux (G) for a 30-min time step ${}^1$. The bold numbers correspond to the best score between both models.

Site and Variable		RMSE ${}^1$ (W m${}^{-2}$)		Pbias (%)		R		NSE
Site	flux	SVS	CLASS	SVS	CLASS	SVS	CLASS	SVS	CLASS
AU-ASM	$LE$	42.0	43.1	3.7	7.1	0.77	0.79	0.43	0.40
	H	62.1	74.2	−17.6	-29.8	0.95	0.96	0.89	0.84
	$R_n$	44.3	48.5	−12.6	-18.2	0.99	0.99	0.97	0.97
	G	38.7	34.9	−157.6	76.9	0.82	0.92	0.66	0.73
US-Wkg	$LE$	71.1	65.7	−10.3	2.6	0.77	0.77	0.28	0.39
	H	61.5	52.4	−5.3	−28.5	0.89	0.91	0.72	0.79
	$R_n$	39.5	46.3	−8.3	−11.8	0.99	0.99	0.98	0.97
	G	56.7	38.2	−204.2	24.0	0.75	0.89	0.52	0.78
Us-Ton	$LE$	40.6	46.6	4.5	16.7	0.83	0.80	0.68	0.58
	H	66.3	79.1	−25.6	−56.5	0.88	0.87	0.74	0.63
	$R_n$	30.0	37.4	−8.8	−20.8	0.99	0.99	0.98	0.97
	G	55.6	68.0	−30.7	55.0	0.65	0.66	−7.01	−10.97
IT-Noe	$LE$	43.2	45.0	0.5	4.8	0.76	0.70	0.17	0.10
	H	64.8	74.9	−15.4	−20.0	0.95	0.95	0.87	0.83
	$R_n$	54.8	47.0	−14.9	−14.4	0.99	0.99	0.96	0.97
	G	43.3	42.8	−100.7	−32.4	0.65	0.82	0.39	0.41
MY-PSO	$LE$	50.2	46.8	−0.8	4.0	0.94	0.95	0.88	0.90
	H	47.9	58.2	−6.7	-24.2	0.89	0.89	0.79	0.69
	$R_n$	38.6	40.3	−5.7	-5.8	0.99	0.99	0.97	0.97
	G	31.4	66.3	−174.2	210.3	0.55	0.35	−62.32	-280.81
GF-Guy	$LE$	75.4	72.1	−12.8	−5.6	0.87	0.88	0.74	0.77
	H	60.4	42.9	86.0	35.1	0.85	0.87	0.21	0.60
	$R_n$	31.4	30.9	−6.1	−6.6	0.99	0.99	0.98	0.98
	G	NA ${}^2$	NA	NA	NA	NA	NA	NA	NA

${}^1$ RMSE, Root Mean Square Error; Pbias, Percent bias; R, Pearson correlation coefficient; NSE, Nash–Sutcliffe efficiency. ${}^2$ NA, Not Available.

.

3.1.1. Arid Sites

The observed annual cycle of $LE$ displays a maximum peak in February during the rainy season at the AU-ASM site (Figure 3a). The $LE$ annual cycle observations are well represented by the SVS simulations, resulting in a good 30-min correlation (R = 0.77, Table 4) and a small overestimation (Pbias = 3.7%, Table 4). H observations peak during the summer months of November–January. Although the mean annual cycle of H is well represented by SVS, it is underestimated by 20 W m${}^{-2}$ on average during the months of October–December. This could be explained by the energy partition at the surface and more precisely by the underestimation of the net radiation (21 W m${}^{-2}$ during the same months) as explained in Section 3.3. $LE$ and H simulations of SVS are similar to that of CLASS, as reflected by the scores in Table 4. $R_n$ is the most dominant term of the annual cycle as well as the best-simulated one in both models (NSE = 0.97, Table 4). However, the bias of $Rn$ is large even if correlation is high. Although incident short and longwave radiations are given, the positive bias suggest that some parameters related to the albedo is unsuitable, as discussed in Section 3.3. In addition, its underestimation during October–December implies that less energy is available for the partition into turbulent fluxes, particularly into H. G is the least dominant flux of the annual cycle, with monthly mean observed value equal to 2.2 W m${}^{-2}$ (Figure 3a). Despite the good correlation of G (R = 0.82 by SVS and R = 0.92 by CLASS, Table 4), SVS shows a large underestimation (Pbias = −157.6%, Table 4) while CLASS shows a large overestimation (Pbias = 76.9%, Table 4). Note that, for the particular case of G, Pbias may lead to large values as a result of the small absolute values of the observations. The conflicting results for G in SVS are associated to: (1) The way G is computed. As the residual of the other energy fluxes, G does not capture properly the seasonality of the observations in the same way than CLASS where G is estimated as a function of soil temperatures. (2) The conceptual fact that the simulations are not really being compared to the point measurement at a flux plate (see Section 3.3 for more details). We obtained similar results for the US-Wkg arid site, which is partially covered by long grass. No valid data are available in January at this site due to the energy closure imposed on the selection criteria.

3.1.2. Mediterranean Sites

Observed $LE$ at US-Ton shows a seasonality with a maximum of 130 W m${}^{-2}$ in May (Figure 3c). The annual cycle of $LE$ is well represented by SVS. SVS produces a value of 70 W m${}^{-2}$ in May and tends to underestimate observations during the growing season (April–June). On the other hand, CLASS tends to better represent the maximum in May (100 W m${}^{-2}$). A reasonable explanation of the poor performance of SVS on May is given by closer inspection of the components of LE, where SVS plant transpiration tends to be lower that CLASS (Figure 4c). H observations show a strong seasonality with a maximum peak of 160 W m${}^{-2}$ in July. SVS simulates this maximum as 155 W m${}^{-2}$. SVS (R = 0.88, Pbias = −25.6%) simulates H better than CLASS (R = 0.87, Pbias = −56.5%, Table 4) at this site. $R_n$ is the best simulated term by SVS (Pbias = -8.8%, NSE = 0.98, Table 4), even better than by CLASS (Pbias = −20.8%, NSE = 0.97). G, with a mean observed value equal to −1.34 W m${}^{-2}$, is poorly simulated by SVS (Pbias = 30.7%, NSE = −7.01) and CLASS (Pbias = −55.0%, NSE = −10.97). The performance of SVS and CLASS is similar at the IT-Noe and the US-Ton site for all energy fluxes.

3.1.3. Tropical Sites

MY-PSO shows a tendency towards uniform $LE$ and H observations throughout the year (see Figure 3e). The behavior of $LE$ (Pbias = −0.8%, R = 0.94, Table 4) and H (Pbias = −6.7%, R = 0.89, Table 4) is well captured by SVS. $R_n$ is equally well simulated, contrary to G (mean observed value equal to −0.87 W m${}^{-2}$), which is largely underestimated (Pbias = 174.2%) and poorly correlated (R = 0.55). SVS simulates all energy fluxes in a similar fashion to CLASS at this site, except for H, which tends to be lower in CLASS (Pbias = −24.2%) than in SVS (Pbias = −6.7%). The $LE$ and H fluxes are accurately simulated by SVS at the GF-Guy site during the wet season (December–August), but are underestimated by SVS during the dry season (September–November), where $LE$ cannot be sustained (Figure 3f). A close examination on the components of $LE$ during the same months (Figure 4e) reveals lower values of SVS canopy evaporation and plant transpiration than CLASS.

3.2. Partitioning of Evapotranspiration

Figure 4 shows the partitioning of evapotranspiration ($ET$) into three different components: soil evaporation ($Es$), transpiration ($Ts$), and wet canopy evaporation ($Ec$). Despite the absence of observations for $ET$, $Es$ and $Ts$, the simulations allow us to describe the variability of each of the component. Note that Figure 4a–f illustrates mean values and not accumulated ones, as is typical of these kind of variables (remember that $ET$ was derived from $LE$, which was filtered according to the energy closure criteria described in Section 2.2).

According to SVS simulations for the AU-ASM arid site, $Es$ is the main $ET$ component ($Es/ET$ = 63%). A similar result was found for the other arid site, US-Wkg, where $Es/ET$ = 84% (Figure 4b). $Ts$ is the main $ET$ component at the US-Ton Mediterranean site ($Ts/ET$ = 63%). The large model differences of maximum $ET$ during the months of April–June are associated with differences in $Ts$ at this site (Figure 4c). The tropical sites illustrate the relevant role of canopy interception and plant transpiration. For example, $Ts/ET$ = 72% at the MY-PSO site, which consists of 70% evergreen forest (Figure 4e). Similarly $Ts/ET$ = 79%, at the GF-Guy tropical site (Figure 4e). The underestimation of $ET$ by SVS during the months of September–November (Figure 4f) is clearly associated with the underestimation of $Ts$, which is probably due to the overestimation of stomatal resistance.

$Ts/ET$ values from the literature can be used to benchmark model results. SVS simulations during the growing season (June–October) at the site US-Wkg lead to a value of 84%, which is above values observed by Moran et al. [50]: values of 45% in 2002 and 75% in 2007. On the other hand, CLASS had a closer estimate of 61% during the same period. A compilation of 81 studies presented by Schlesinger and Jasechko [84] partitioned annual $ET$ into $Ts$ at the ecosystem scale in Mediterranean (47 ± 10%) and tropical ecosystems (70 ± 14%), which are in line with SVS and CLASS simulations at the Mediterranean (respectively, 50% and 57% on average) and tropical (respectively, 76% and 60% on average) sites.

3.3. Sources of Uncertainties Associated with Energy Fluxes

SVS surface energy flux simulations show a tendency to underestimate $R_n$ at all sites. More precisely, the average underestimation in terms of the Pbias score is of the order of 10.4% at the arid sites (AU-ASM, US-Wkg), 11.9% at the Mediterranean sites (US-TON, IT-Noe) and 5.9% at the tropical sites (MY-PSO, GF-Guy). Alves et al. [30] showed the same scale of underestimations in $R_n$ when they contrasted CLASS simulations with five FLUXNET sites representing semi-arid, temperate, and tropical climate types. The main reason is the underestimation of net shortwave radiation for which the percent bias score is on the same order as $R_n$: 10.3% at arid sites, 13.3% at Mediterranean sites and 3.5% at tropical sites. Furthermore, albedo simulated values in SVS (and also CLASS) tend to be overestimated at all the studied sites. Indeed, daily simulated albedo, averaged between 10:00 and 14:00, for SVS(CLASS) showed poor performances: RMSE = 0.07, Pbias = 60.28, r = 0.24 (RMSE = 0.07, Pbias = 48.23, r = 0.16). Part of the bias may be related to instrumentation uncertainty in the FLUXNET database. Wilson et al. [85] showed that $R_n$ observations in FLUXNET are typically characterized by systematic instrument errors of up to 15%.

For turbulent fluxes, H at arid (AU-ASM, US-Wkg) and Mediterranean sites (US-Ton, IT-Noe) is better simulated by SVS (NSE = 0.80 on average) than $LE$ (NSE = 0.39 on average), while the opposite is true at the tropical sites (MY-PSO, GF-Guy) where $LE$ (NSE = 0.81 on average) is better simulated than (NSE = 0.50 on average). The fact that $LE$ is less performing than H at arid and Mediterranean sites is because of the presence of larger variability in biophysical characteristics, where significant climate seasonality exists. A closer look at the evapotranspiration components (Figure 4) reveals that transpiration is the most likely process to be biased (here it is essential to mention the benchmarking role of the CLASS model). Another factor influencing the simulation of turbulent fluxes may be the nocturnal observations. At nearly all FLUXNET sites, it is assumed that the advection of scalars can be neglected [86]. However, non-zero values of mean vertical velocity and vertical advection are possible and can be induced by horizontal heterogeneity or mesoscale and synoptic scale forcings [87,88].

G simulated by SVS is biased for all sites without exception; one reason may be the method associated with the estimation of soil temperature. As reported by Verseghy [14], the force-restore method causes the surface temperature to respond to diurnal forcing rather slowly, leading to a damping of G (and of turbulent fluxes). Additional studies concerning the implementation of force-restore method into LSMs showed accurate estimates of shallow (upper) soil layers (e.g., [40,89,90,91,92]), but limitations on deep soil temperature, and under snow cover and frozen soils [93,94]. Furthermore, soil heat flux measurements using plates, which is typical in FLUXNET, can be inaccurate because the thermal conductivity of the heat flux plate and the surrounding soil may be unequal [95]. Heat flux plates typically have a conductivity similar to that of soil, but thermal properties vary temporally because of changes in soil water content and temperature [96]. Soil heat flux plates can also alter the environment they are measuring, notably by limiting water movement in soils [97].

3.4. Water Balance

Figure 5 reveals the contribution of the different water balance terms from Equation (2) for the various studied climate conditions when precipitation is considered to be 100%. According to the SVS simulations, $ET$ dominates with 93%, 84% and 53% in arid, Mediterranean and tropical sites, respectively (Figure 5a–f). The remaining contribution is split mainly between overland ($R_{surf}$) and base ($R_{base}$) flows. Although CLASS simulates $ET$ similarly to SVS, differences between the distribution of $R_{surf}$ and $R_{base}$ are observed. For example, at the MY-PSO tropical site, $ET$ is around 65% of precipitation in both models. However, $R_{surf}$ is 35% in SVS, but only 8% in CLASS (Figure 5e). Such a difference between models could lead to a substantial difference in terms of hydrological responses if $R_{surf}$ and $R_{base}$ are not aggregated together in the routing model. This is the case in GEM-Hydro, where $R_{base}$ is first provided to an intermediate reservoir in the routing scheme Watroute (aiming to represent the aquifer), before reaching the channel [6].

CLASS guarantees water balance closure for each of its four subareas up to the accuracy limit of 0.1 kg m${}^{-2}$ [98], as this is a requirement of long-term simulations of climate models [72]. SVS closes the water balance with a small positive bias (values in parenthesis) at all sites: AU-ASM (0.2%), US-Wkg (0.2%), US-Ton (0.6%), It-Noe (0.1%), MY-PSO(0.05%) and GF-Guy(0.04%).

3.5. Soil Moisture

Scores obtained from the evaluation of the soil water content simulations are summarized in

Table 5. Scores between observations and simulations (SVS and CLASS) of the surface water content for all the study sites. The scores were obtained for 30-min sampling intervals. The bold numbers correspond to the best score.

Site	RMSE ${}^1$ ( m${}^3$ m${}^{-3}$)		Pbias (%)		R		NSE
Site	SVS	CLASS	SVS	CLASS	SVS	CLASS	SVS	CLASS
AU-ASM	0.03	0.03	35.1	40.2	0.92	0.93	−0.04	0.20
US-Wkg	0.03	0.02	−14.2	−8.0	0.87	0.91	0.64	0.78
US-Ton	0.07	0.09	−22.9	−30.4	0.93	0.92	0.73	0.60
IT-Noe	0.08	0.10	−23.3	−30.2	0.84	0.83	0.11	−0.38
MY-PSO	0.17	0.15	−39.4	−36.2	0.86	0.86	−26.50	−22.71
GF-Guy	0.14	0.12	72.1	57.1	0.78	0.76	−5.60	−3.74

${}^1$ RMSE, Root Mean Square Error; Pbias, Percent bias; R, Pearson correlation coefficient; NSE, Nash–Sutcliffe efficiency.

. They were computed by comparing 30-min point-scale measurements of the shallowest layer (typically measured at a 5- or 10-cm depth) to modeled values corresponding to the topsoil layer (0–10 cm). For deeper layers (10–410 cm), only SVS and CLASS simulations are available and are presented in Figure 6, Figure 7 and Figure 8.

3.5.1. Arid Sites

The surface layer soil moisture content is highly modulated by precipitation at all arid sites. At AU-ASM (Figure 6a), SVS surface soil moisture performance is low (NSE = −0.04, Table 5) and less than CLASS (NSE = 0.20, Table 5). During the dry periods (from May to September), SVS tends to represent the minimum observed value of 0.05 m${}^3$ m${}^{-3}$ properly. Exceptions occurred during the dry periods of 2013 and 2014 when SVS depleted soil moisture to 0.0 m${}^3$m${}^{-3}$ while CLASS was comparable to observations (0.05 m${}^3$ m${}^{-3}$). As reported by Maheu et al. [46], SVS depletion during dry periods at semiarid sites is explained by the absence of the residual soil water content, which in CLASS was defined as 0.04 m${}^3$m${}^{-3}$ [14]. During rainy events, SVS overshot the observed maximum soil water content (Pbias = 35.1%, Table 5). For deeper layers (10–420 cm), SVS is slightly wetter than CLASS. As a result, the minimum soil water content of 0 m${}^3$ m${}^{-3}$ is reached at the 3-m depth in SVS, while it occurred at the 2.4-m depth with CLASS (Figure 6c,d). As a consequence, more water is available in the soil and is transferred as overland flow, drainage, or as plant uptake through the roots. Figure 5 shows around 5% more overland flow for SVS than for CLASS at this particular site. A good representation of near-surface soil moisture by SVS and CLASS (NSE = 0.64 and 0.78, respectively, Table 5) is observed for the other arid site US-Wkg (not shown).

3.5.2. Mediterranean Sites

SVS performed well at the US-Ton site during the study period (NSE = 0.73, Table 5), even better than CLASS (NSE = 0.60, Table 5). During the wet season, SVS systematically underestimated surface water content, the maximum observed soil water content was 0.55 m${}^3$m${}^{-3}$, while the maximum simulated was 0.45 m${}^3$m${}^{-3}$. SVS also slightly underestimated surface soil moisture during the dry seasons. This response is reflected by a negative Pbias (=−22.9%, Table 5). Soil water content for SVS shows significant water percolation (soil moisture between 0.3–0.4 m${}^3$m${}^{-3}$) during the wet season and a progressive water transfer from the surface layer during the dry season (see Figure 7c). This result is contrasted to that of CLASS, which shows slow water percolation during the wet season and drastic transitions during the dry season, resulting in generally drier conditions (Figure 7c). Considering that both models account with the same meteorological forcing and soil texture, divergences in the deep soil moisture could be explained by differences in the parameter setting as a function of the soil texture (Table 2). See also Section 3.6 for a more detailed discussion.

3.5.3. Tropical

Both models showed low performance for surface water content simulations at the MY-PSO tropical site (NSE = −26.50 and −22.71 for SVS and CLASS, respectively, Table 5). SVS underestimates water content mainly during the dry period. This behavior is depicted by Pbias = −39.4%. Despite the weak performance of SVS, its similarity to CLASS and the reasonably good correlation of both models (R = 0.86 for SVS and CLASS, Table 5) suggest that there is some bias associated with the point-scale measurements that may be associated with local inaccuracies in the soil texture. In fact, the two primary studies associated with the site [54,99] represent soil moisture as the average at 0.1, 0.2 and 0.3 m with a closer variability to our simulations. Additionally, there was large variability in the datasets employed by Cosby et al. [80] when the relationships between soil texture and hydraulic parameters were developed. For deeper layers, both models tend to simulate uniform wet conditions (see Figure 8c,d), which is contrary to the arid and Mediterranean sites.

3.6. Uncertainties Associated with Soil Moisture

A mismatch can occur when we compare point-scale soil moisture observations against the averaged simulated values of the top 10 cm layer, which is somewhat close to the surface. Furthermore, the fact that the soil texture had to be defined as uniform over the entire soil column, as already mentioned, certainly has an additional influence on the representativeness of the results. On tropical sites with evergreen tropical forest, the soil surface is usually shaded and covered with organic-matter litter that is frequently wet [100]. It imposes a condition that reduces the surface runoff and increases the base flow [100], which may explain the divergences between SVS and CLASS water balances at tropical sites (Figure 5e,f). Concerning model structure, SVS and CLASS models are not directly informed on the local soil hydraulic properties, which incidentally are unknown at the selected sites. Soil hydraulic properties are instead deduced from pedotransfer functions that provide the most probable values for a given soil texture (Table 2). This procedure inevitably leads to some divergence with reality that will push the models to diverge from observations. In addition, although both models are similarly process-based as to soil description, the constant values differ in the parameterized equations. This is the case for the saturated hydraulic conductivity, which provides higher values in SVS than in CLASS, allowing a greater water content in deeper layers as well as a role in the water uptake from the soil. This effect is observed at the US-TON site (Figure 5). At the same time, model discrepancies in the description of stomatal resistance may also have a role in water uptake from the soil. A closer introspection of the model structure bias may require a sensitivity analysis of specific parameters, particularly the hydraulic conductivity and stomatal resistance.

4. Conclusions

This study used flux tower observations to assess the ability of the SVS model to simulate surface energy fluxes and surface soil moisture under snow-free conditions for various climates and biomes. In addition to observations, the CLASS model was used as a benchmark. This approach provided an integral diagnostic of SVS in terms of energy and water exchange between the land surface and the atmosphere.

SVS led to realistic simulations of latent heat fluxes (NSE = 0.58 on average), sensible heat fluxes (NSE = 0.70 on average) and net radiation (NSE = 0.97 on average) under arid (sites AU-ASM and US-Wkg), Mediterranean (sites US-Ton and IT-Noe) and tropical (sites MY-PSO and GF-Guy) climates. Soil heat fluxes were well simulated at the AU-ASM (NSE = 0.66) and US-Wkg (NSE = 0.52) arid sites, reasonably well simulated at one of the Mediterranean sites (IT-Noe, NSE = 0.39) and poorly at one of the tropical sites (MY-PSO, NSE = −62.32).

Surface soil moisture was well simulated at one arid (US-Wkg, NSE = 0.64) and one Mediterranean site (US-Ton, NSE = 0.73), poorly simulated at one arid (AU-ASM, NSE = −0.04) and one Mediterranean site (IT-Noe, NSE = 0.11) and very poorly simulated at the tropical MY-PSO (NSE = −26.50) and GF-Guy (NSE = −5.60) sites. The simulations of SVS and CLASS showed a consensus for deep soil water content at arid and tropical sites, but SVS showed significant water percolation during the wet season and slow percolation during the dry season for Mediterranean sites, while CLASS produced slow water percolation during the wet season and drastic and dry transitions during the dry season.

Evapotranspiration was the most significant term for water balance with contributions of 95%, 82% and 50% in arid, Mediterranean, and tropical sites, respectively.

As a result of this study, some deficiencies in the SVS model’s ability to simulate energy fluxes and surface soil moisture were identified and thus deserve further attention:

• The albedo of vegetation is simply based on a weighted mean of look-up table values of fixed albedos. This scheme may lead to a negative bias of up to 12% in the net shortwave radiation and 13% in the net radiation. Albedo is specified the same way as root depth, density and LAI. These are all parameters that can be improved and will impact SVS simulations. The photosynthesis also has several look up tables for various key parameters controlling the stomatal resistance. External databases and satellite-derived approximation could help better specify some of these parameters.
• Abrupt underestimations of $LE$ at GF-Guy site is mainly associated with transpiration simulations. Thus, further attention should be given to this process.
• Poorly simulated G is associated with the simple method used (residual of the other energy fluxes) and the single-layer description of the soil. To avoid these structural biases, vertical transport of heat in the soil along with multilayer-energy balance would be an asset in the soil description.
• The better simulation of soil moisture with CLASS in arid sites suggests that certain features, such as the inclusion of residual saturation when defining water retention curves, can be advantageous. Alternatively, limits on evaporation and moisture stress parameterizations may help to prevent soil moisture to reach zero values.
• Surface soil moisture at tropical sites should take into account organic matter.