From Standard Weather Stations to Virtual Micro-Meteorological Towers in Ungauged Sites: Modeling Tool for Surface Energy Fluxes, Evapotranspiration, Soil Temperature, and Soil Moisture Estimations

Abstract

One of the benefits of training a process-based, land surface model is the capacity to use it in ungauged sites as a complement to standard weather stations for predicting energy fluxes, evapotranspiration, and surface and root-zone soil temperature and moisture. In this study, dynamic (i.e., time-evolving) vegetation parameters were derived from remotely sensed Moderate Resolution Imaging Spectroradiometer (MODIS) imagery and coupled with a physics-based land surface model (tin-based Real-time Integrated Basin Simulator (tRIBS)) at four eddy covariance (EC) sites in south-central U.S. to test the predictability of micro-meteorological, soil-related, and energy flux-related variables. One cropland and one grassland EC site in northern Oklahoma, USA, were used to tune the model with respect to energy fluxes, soil temperature, and moisture. Calibrated model parameters, mostly related to the soil, were then transferred to two other EC sites in Oklahoma with similar soil and vegetation types. New dynamic vegetation parameter time series were updated according to MODIS imagery at each site. Overall, the tRIBS model captured both seasonal and diurnal cycles of the energy partitioning and soil temperatures across all four stations, as indicated by the model assessment metrics, although large uncertainties appeared in the prediction of ground heat flux, surface, and root-zone soil moisture at some stations. The transferability of previously calibrated model parameters and the use of MODIS to derive dynamic vegetation parameters enabled rapid yet reasonable predictions. The model was proven to be a convenient complement to standard weather stations particularly for sites where eddy covariance or similar equipment is not available.

1. Introduction and Goals

The magnitude of the surface energy fluxes including net radiation (NR), latent heat flux (LE), sensible heat flux (H), and ground heat flux (G) impact both atmospheric and land surface processes in relation to water, energy, and biogeochemical cycles [1,2,3,4]. Furthermore, evapotranspiration, soil temperature, and moisture are key in triggering and maintaining drought and flooding conditions. In all these processes, vegetation cover, type, and activity are key factors that control, among others, transpiration rates, ground albedo, below-canopy soil temperature and moisture, and solar radiation sheltering [5,6,7,8]. Annual phenologic changes in vegetation activity and cover result in shifts in partitioning of the surface energy balance (SEB), alterations in the surface and root-zone soil water content, and changes in the land surface temperatures that are key in the land–atmosphere interactions [9,10,11,12]. Furthermore, sudden natural or human-induced changes in land cover, vegetation type, and canopy biomass influence the SEB fluxes, including abrupt shifts in evapotranspiration (ET), soil moisture, and ground temperature that then influence regional hydro-climatic patterns [13,14] and long-term water, carbon, and energy fluxes and stocks [15,16,17,18]. Therefore, the inclusion of dynamically changing vegetation parameters like albedo, stomatal resistance, and leaf area index, among others, can significantly improve the accuracy of any land surface model.

The eddy covariance (or eddy correlation, EC) method was proposed from measurements of mass flux and momentum to compute energy, vapor, and carbon fluxes within the atmospheric boundary layer [19,20,21]. The theory supporting the EC method poses that, under conditions of horizontal homogeneity, the net transport between the surface and the atmosphere exists in one dimension. Due to this, the method can be applied to estimate the flux density between turbulent fluctuations of the variable of interest and the vertical wind [22]. Besides some limitations regarding the accuracy of some energy exchange sensors [22], the EC method continues to be widely used for energy, vapor, and carbon flux studies around the world [23]. Usually, energy, water, and carbon fluxes are accompanied by measurements of the soil water content and temperature at different depths and standard weather variables. However, despite the multiple benefits for science and engineering that these EC systems provide as micro-meteorological arrays, the cost and complexity of purchasing, maintaining, and collecting the data constrain our understanding of the physical processes that rule various ecosystem types across different latitudes [24]. Due to this limitation, although not comparable to the appropriateness of an EC system, some alternative techniques have been used to quantify SEB components, many of them focusing exclusively on latent heat flux and thus evapotranspiration (ET).

For estimation of the SEB and ET, some statistical and empirical approaches based on remote sensing imagery have been proposed [25] including methods applicable over large areas [26]. Some models such as the Surface Energy Budget Algorithm for Land (SEBAL [27,28]), the Simplified Surface Energy Balance Index (S-SEBI [29]), the Mapping ET with Internalized Calibration (METRIC, [30]), the Surface Energy Balance System (SEBS [31]), and the operational Simplified Surface Energy Balance (SSEBop [32]) are promising but often constrained to the spatiotemporal resolution of the remotely sensed imagery and lack the inclusion of micro-meteorological and terrain factors. Furthermore, besides their complexity, the results tend to show overestimations of sensible, ground, and latent heat flux values during dry periods or underestimations of ET during wet conditions [25,33], suggesting the necessity to better incorporate dynamic vegetation processes [25]. Another approach uses the MOD16 model global evapotranspiration project for which the improved algorithm [34] estimates LE (and thus ET) by using the Pennman–Monteith equation for a reference crop. However, besides the temporal (eight-day) and spatial resolution (250 m) limitations, there are important biophysical parameters, such as the mean potential stomatal conductance, that despite changing across biomes [35], are assumed to be constant in time [34]. This assumption becomes a model limitation, for example, with plant species that use the crassulacean acid metabolism (CAM) and open their stomata at night to reduce water loss (common in hot and dry areas), therefore contributing to uncertainties in the ET response up to 25% at the catchment scale [36]. Physics-based approaches often work at large spatial scales with low pixel resolution (e.g., WRF-Hydro at 1 km) that end up aggregating important spatial variability into the final coarse-pixel results. Furthermore, land surface models often simplify the role of vegetation to a static component, prescribing characteristics based on topography or climate [37,38]. Such approaches often over-parameterize vegetation-related processes. For example, the inclusion of dynamic vegetation requires additional multi-parameterization of the land surface model Noah-MP [39] embedded within WRF-Hydro [40], making the modeling process more complex and less efficient.

1.1. Fusing Remote Sensing with a Multi-Physics Framework: The Key to Model Transferability

To date, the three grand challenges of land surface models are (1) managing process complexity, (2) representing land surface heterogeneity, and (3) understanding parametric dynamics [41,42]. Together, they introduce significant uncertainty and challenges in tractably representing processes. The explicit consideration of vegetation dynamics into water and energy budgets at the critical zone tackles the third challenge as an unavoidable condition for continuous-in-time and robust model simulations of land–atmosphere interactions [41,42]. Nonetheless, informing vegetation model parameterization from field observations is challenging and made difficult by the lack of in situ data about plant’s phenology, activity, and cover. This is one of the reasons why, to date, hyper-resolution models end up avoiding vegetation dynamics and rather opt for extensive calibration as an inverse (and often leading to convoluted equifinality) approach to intentionally improve simulation skill scores. Despite the boom in the use of remotely sensed information, linkages between spectral values and the parameters that drive energy and water fluxes in vegetation are seldom connected with remote sensing through multi-physics frameworks all at once. This article presents and tests a framework for vegetation processes in relation to shortwave and longwave energy fluxes, precipitation interception, vegetation shading, and transpiration and their relations to remotely sensed information from Moderate Resolution Imaging Spectroradiometer (MODIS) (MCD15A3H, MCD43A4, MCD15A3H, and MCD43A) for sub-weekly and daily updates on leaf area index (LAI), Normalized Difference Vegetation Index (NDVI), photosynthetically active radiation (PAR), and albedo. These important parameters control key reservoirs and fluxes of water and energy such as the maximum vegetation canopy field storage capacity (S), water throughfall coefficient (p), vegetation optical transmission coefficient (k${}_t$), vegetation stomatal resistance (r${}_s$), vegetation cover (v${}_f$), and the amount of absorbed photosynthetically active radiation by plants (Q). The ability to estimate each of those parameters accurately and to provide their time series facilitates not only the performance of any land-surface model but also improves the physical meaning, understanding, and tractability of internal processes. In places where soils and vegetation may be similar to EC sites at which the model is trained, the transferability of such a model would be conveniently facilitated by the use of these parameter-inference equations that rely only on remote sensing and use previously calibrated (from similar sites) static model parameters.

1.2. Goals

This study seeks to evaluate the capability of a process-based land surface model working at the eddy covariance footprint scale to estimate net radiation and SEB fluxes along with soil moisture and temperature and to test the model transferability to hydrologically similar sites without parameter calibration or recalibration. The model was run over four EC sites in Oklahoma, paired according to vegetation type (i.e., grassland and cropland) considering dynamic vegetation changes in a computationally efficient format. Section 2 presents the materials and methods, while Section 3 summarizes the results including model calibration, validation, and framework transferability to similar sites. A scientific discussion is provided in Section 4, and the conclusions are synthesized in Section 5.

2. Data and Methods

2.1. Study Sites

The process-based model was calibrated, validated, and tested at four EC system sites within the north-central plains of the state of Oklahoma (see Figure 1).

Due to their geographic proximity, the climate and weather patterns are similar among the four sites. Within the region, during the study periods, the annual average temperature and precipitation values were 15.1 ${}^{\circ }$C and 503.2 mm/y, respectively [43]. Oklahoma has more than 492 vegetative species [44] with the most common being tree, shrub, and herbaceous vegetation types. Cross timbers, central mixed grass, high plains short grass, and Osage tall grass are the most common grassland types, covering over a third of the state area [44]. Cropland, and deciduous and evergreen forests cover most of the state area (Figure 1). The remainder of the state land is covered by species that make up less than 10,000 ha each [44]. The two selected vegetation types for this study (i.e., grassland and cropland) are representative of the dominant land cover patterns of the rural grounds of the north-central section of the state. Within a typical year, vegetation cover and activity respond to phenologic changes, and Phenocam photos taken at the study sites in contrasting cool and warm-season months clearly illustrate such phenologic changes (see Figure 2).

Table 1. Geographic coordinates, soil and vegetation type, and purpose of the four selected EC sites.

ID	Lat, Lon	Soil and Vegetation Type	Purpose Within the Study
ARM-CF	36.6058N, 97.4888W	Silty clay loam. Crop field (winter wheat, soy, corn, and alfalfa)	Model calibration and validation in cropland
MOISST	36.0634N, 97.2169W	Sandy clay loam. Rangeland with grazed cattle pasture	Model calibration and validation in grassland
ARM-A74	36.8084N, 97.5488W	Silt Loam. Croplands and rotational crops (i.e., soybean and corn) followed by harvest and a bare soil period	ARM-CF parameter transferability evaluation in cropland
ARM-A32	36.8192N, 97.8197W	Kirkland silt loam. Grasslands (Medford hay pasture) periodically cut for hay	MOISST parameter transferability in grassland

describes the geographic coordinates, soil and vegetation type, and purpose of inclusion (in light of the objectives of this study) of each study site. The sites were paired by vegetation type to ensure similarity for model parameter transferability and testing despite there being some differences between the selected pairs, particularly with regard to management (e.g., prescribed fire during some years at the MOISST site but not at US-A32).

The EC systems shown in Figure 1, excluding MOISST, have their data available under the AmeriFlux network website that provides free and quality-controlled 30 min soil, micrometeorological, energy, water, and carbon flux data. The ARM-CF site serves as a ground validation site of MODIS_L3 products, and the MOISST site is managed by Oklahoma State University (OSU) and is a benchmarking location for the evaluation of in situ soil moisture sensing technologies. At MOISST, sections of the field are burned for woody plant control every three years.

The areal extent of the model domain was determined by the size of the eddy flux source areas (i.e., fetch or footprint) driven by the height of the flux sensors (2 m above ground), topographic and vegetation roughness, and aerodynamic conditions. The numerical computation of the eddy flux footprint areas (Figure 3) was conducted using a Lagrangian model for various turbulent stratifications with backward paths [45] that required a 2D parameterization across sites. Figure 3 shows that flux footprints are roughly elliptical, with a south–north major axis across stations. Nonetheless, the diameter of 80% of the contributing area is no wider than 50 m from the center tower in all cases.

2.2. Data

2.2.1. Terrain and Vegetation

Topography, soil, and vegetation properties determine the surface water and energy interactions, the soil control volume and its properties, and the model parameter and boundary conditions. The flux footprint of each EC system (Figure 3) was used to determine the maximum extent of terrain and vegetation input files.

Table 2. Topography, soil, and vegetation data sources with spatiotemporal resolution.

Type of Input	Source	Product	Spatial Resolution	Temporal Resolution
Digital Elevation Model	USGS	SRTM	30 m	N.A.
Soil type	Logs, NRCS	Texture	Footprint	N.A.
Land cover type	USGS	NLCD	30 m	2016 version
Leaf Area Index (LAI)	MODIS	MCD15A3H	500 m	4 days
NDVI	MODIS	MCD43A4	500 m	daily
Photosynthetically Active Radiation (PAR)	MODIS	MCD15A3H	500 m	4 days
Albedo	MODIS	MCD43A	500 m	daily

compiles model input type, source, and spatial and temporal resolutions. A digital elevation model was extracted from the USGS Shuttle Radar Topography Mission (SRTM) for the simulation footprints (see Figure 3). Soil textural types were obtained from each site’s metadata, from site surveys, and from the Natural Resource Conservation Service (NRCS). Vegetation types were obtained from the USGS National Land Cover Dataset (NLCD). A number of initial (precalibration) static parameters for soil and vegetation were retrieved from this information (see Section 2.4 for parameter value tuning). However, dynamic vegetation parameters were obtained from the time series of Leaf Area Index (LAI), Normalized Difference Vegetation Index (NDVI), and surface albedo (Al) values from MODIS since it provided the best combination of spatial, temporal, and spectral resolutions.

2.2.2. Weather Forcing

Hourly weather forcing data were used as input to the model. Incoming shortwave radiation (SW), air temperature (T), air vapor pressure (VP), wind speed (WS), atmospheric pressure (Pa), and precipitation (P) forced the model at each study site. A careful description of the type of sensors measuring each of these variables at EC sites can be found at the FLUXNET website and [46]. Since weather forcing must be continuous in time, a few missing data were interpolated to guarantee completeness.

2.2.3. Surface Energy Fluxes, Surface, and Root-Zone Soil Temperature and Moisture

The hourly data that were used to calibrate and evaluate the land surface model, with regard to the SEB fluxes, were net radiation (NR), latent heat flux (LE), sensible heat flux (H), and ground heat flux (G). The EC system instrumentation is standard and consists of a RAD-Net radiometer for NR, a RAD-SW Pyranometer Class2 for the incoming short wave radiation, a GA OP-LI-COR LI-7500 to measure LE, a SA-Gill Windmaster Pro for H, and ground heat plates for G [46]. The measured fluxes were not corrected for energy balance closure as physics-based simulation equations do not rely on closing the balance either, thus avoiding the introduction of biases [47,48]. Additionally, the soil-related predictors were surface soil moisture (SSM), root-zone soil moisture (RSM), surface soil temperature (SST), and root-zone soil temperature (RST). SSM, SST, RSM, and RST were measured using Stevens Water Hydra Probe II sensors [49]. The measurements at the root-zone level were taken at approximately 0.9 m below the ground surface at the MOISST and at 0.3 m at US-A32 and US-A74. A detailed description of the sensors used to measure soil-related variables can be found in [50].

2.3. Modeling Framework

2.3.1. Multi-Physics Model

tRIBS is a process-based, distributed hydrologic model developed at the Ralph M. Parsons Laboratory, Massachusetts Institute of Technology [47,48,51,52]. It uses time evolving, spatially distributed atmospheric, land cover, soil, and topographic information at each Voronoi element to simulate both energy and water flows within a computational domain. Elements are connected through surface and subsurface water and energy fluxes throughout the domain. The process-based framework includes parameterizations of the surface and subsurface processes such as net radiation, heat fluxes, rainfall interception, evapotranspiration, infiltration with continuous soil moisture and temperature accounting, lateral moisture transfer in the unsaturated and saturated zones, and hillslope and channel runoff routing [47,48].

Table 3. Tin-based Real-time Integrated Basin Simulator (tRIBS) model of the physics of energy fluxes and soil status including mathematical framework and reference sources.

Symbol	Description	Method	Reference
NR	Net Radiation	Based on the four vertical components of the radiation budget at the surface including incoming and outgoing short- and longwave components NR = R${}_{si}$ + R${}_{li}$ − R${}_{so}$ − R${}_{lo}$. All terms are computed from standard weather (e.g., T and VP), surface (SST), and remote sensing measurements (albedo and LAI).	[62,63,64]
LE or ET	Latent Heat Flux or Actual Evapotranspiration	Using the Penmann–Monteith approach, the model partitions reference ET among evaporation from soil, and evaporation from vegetation interception and transpiration. Estimated actual ET accounts for soil moisture as a limiting factor when atmospheric demand is high; wind speed, water vapor deficit, vegetation height, vegetation cover (from LAI), and activity (from NDVI) that determine optical transmission; and atmospheric and stomatal resistances.	[65,66,67,68,69]
H	Sensible Heat Flux	Uses an aerodynamic resistance approach between surface and air temperatures. The atmospheric resistance term depends on wind speed and rugosity terms.	[70]
G	Ground Heat Flux	Based on a force-restore method that solves the heat diffusion equation between soil surface and deeper layers. The flux G is obtained from G = 0.5·C${}_s$d${}_1$($\xi$(d(SST)/dt) + $\omega$(SST-RST)), where C${}_s$ is the soil heat capacity, $\omega$ is the daily frequency of oscillation, d${}_1$ = (2k/$\omega$) is the soil heat damping depth, k = k${}_s$/C${}_s$ is the soil diffusivity, and k${}_s$ is the soil heat conductivity (see Table 5). $\xi$ is computed using Hu and Islam (1995) parameterization.	[71,72]
SSM and RSM	Surface and root-zone soil moisture	A ponding and infiltration scheme based on the kinematic approximation for unsaturated flow for a sloping, heterogeneous anisotropic soil. A soil moisture state results from infiltration, runoff, and subsurface flows and is coupled to loses from soil evaporation and transpiration. The model considers ponded infiltration, infiltration under unsaturated conditions, wetted wedge dynamics for the unsaturated phase, and perched zones and keeps track of the evolution of fronts. Surface and root-zone moisture are integrated within the first 5 cm of soil and at 1 m depth. Soil water content is expressed as a fraction of the soil porosity or degree of saturation.	[73,74,75]
SST and RST	Surface and root-zone soil temperature	SST and RST are obtained during calculation of the transient-state energy budget equation at the surface, C${}_s$·(d(SST)/dt) = Rn-LE-H-G, and the calculation of G (see above). Soil heat wave damping depth and damping depth temperature are intrinsically computed when resolving with the force-restore method to calculate G.	[12,71,72]

summarizes the physics-based equations and literature references that explain the mathematical approaches of each of the variables that are simulated in this study using tRIBS. tRIBS has previously been used for energy, carbon, and water budget simulations at the plot [53,54] and watershed scales [17,40,55,56,57,58,59,60,61] including soil moisture and temperatures. For a detailed description of the equations and physical framework, please refer to Ivanov et al. (2004). tRIBS is an appropriate model for the type of science question tackled in this study due to its high spatial resolution and computational efficiency.

2.3.2. Model Training and Verification

tRIBS was calibrated and validated at the ARM-CF and MOISST sites during periods of time involving cold and warm seasons. The selection of these time intervals was based on data availability and a low number of information gaps.

Table 4. Selected time periods for tRIBS model calibration, validation, and parameter transfer assessment.

Station	Purpose	Simulated Interval	Simulated Period (h)
ARM-CF	Calibration	04/30/2004–06/29/2005	10,000
ARM-CF	Validation	07/01/2008–07/01/2009	8780
MOISST	Calibration	12/09/2013–09/10/2014	7300
MOISST	Validation	11/09/2015–29/12/2016	10,000
ARM-A74	Transferability Evaluation	01/01/2016–06/01/2017	12,384
ARM-A32	Transferability Evaluation	01/01/2016–06/01/2017	12,384

illustrates the time intervals for each of the model training and evaluation procedures at each of the sites. In all cases, a spin-up time of 200 h was allowed to obtain a close-to-steady-state condition for the underground water storage during an initial dry (no rain) period.

Model calibration was conducted for the static soil and vegetation parameters described in

Table 5. tRIBS static calibrated soil and vegetation (underlined) parameters.

Parameter	Description	Parameter	Description
*K${}_S$ (mm/h)	Saturated hydraulic conductivity	${\theta }_s$ (unitless)	Soil Moisture at saturation
${\theta }_r$ (-)	Residual soil moisture	m (-)	Pore distribution index
${\mathsf{\Psi }}_b$ (mm)	Air-entry pressure	f (unitless)	Conductivity decay with depth
A${}_s$ (-)	Saturated anisotropy ratio	A${}_u$ (-)	Unsaturated anisotropy ratio
n (-)	Soil porosity	k${}_s$ (J/msK)	Soil volumetric heat conductivity
C${}_s$ (J/m${}^3$K)	Soil heat capacity	K (mm/h)	Vegetation throughfall drainage coefficient-Rutter
b${}_2$ (mm${}^{-1}$)	Vegetation throughfall drainage exponential parameter-Rutter	H (m)	Vegetation height
${\theta }^{*}$ (-)	Evaporation stress threshold for soil evaporation (-)	${\theta }_s$ (-)	Stress threshold for plant transpiration

*K${}_s$ is at the soil surface and decreases with depth according to f.

, while dynamic vegetation parameters were updated according to the remote sensing information and physics-based equations shown in

Table 6. Physics-based equations linking remote sensing with tRIBS dynamic vegetation parameters.

Parameter	Equation	Remarks
Canopy Field Capacity-Rutter (S, mm)	S = 0.5.LAI	Controls depth of rainfall interception as a function of LAI [79]. The values can range among ecosystems (i.e., SGP OK 0.8–1.2 mm) [47].
Free throughfall coefficient-Rutter (p, -)	p = e${}^{(-1.5LAI)}$	Drives the fraction of rainfall not captured by plants as a function of LAI [54,79].
Optical transmission coefficient (k${}_t$, -)	k${}_t$ = e${}^{-k.LAI}$	Based on Beer–Lambert law. k is the light extinction coefficient determined from [80].
Minimum stomatal resistance (r${}_s$, s/m)	r${}_s$ = $\frac{\frac{Q_{50}+Q}{g_{max}.Q}}{LAI}$	Based on the energy-limited relation by [35,81]. Q${}_{50}$ is the value of the absorbed photosynthetically active radiation (Q) when the maximum seasonal stomatal conductance (g${}_{max}$) is half of its value. LAI is used to upscale the individual leaf estimation to the entire canopy [82].
Absorbed photosynthetically active radiation (Q, W/m${}^2$)	Q = 0.45 SW fPAR	Q drives photosynthesis and transpiration. fPAR is the fraction of photosynthetically active radiation absorbed by plants; 0.45 is the fraction of shortwave (SW) radiation used during photosynthesis [83].
Vegetation Fraction (v${}_f$, -)	v${}_f$ = $\frac{NDVI-NDV{I}_{min}}{{(NDV{I}_{max}-NDV{I}_{min})}^2}$	Vegetation fraction computed as a function of NDVI based on [84]. v${}_f$ plays a determinant role in model transpiration [54,85].
Albedo ($\alpha$, -)	$\alpha$	Absolute value of ground reflectivity.

. Without exception, all parameters shown in Table 5 were calibrated using the procedure explained next. First, a One-At-a Time (OAT) sensitivity analysis was conducted to determine the most important static parameters for each simulation site (Table 5) [17,57]. After this, the Shuffled Complex Evolution (SCE) algorithm was applied to achieve the automated calibration of selected soil parameters [76]. Comparisons of observed and simulated values of NR, LE, H, G, and SST were used to calibrate the model at ARM-CF. At MOISST, NR, LE, H, SST, SSM, and RSM were the variables used to conduct the model calibration. These selections were made based on observed time-series quality, or the lack thereof, at each site. The calibration objective functions were based on the Pearson correlation coefficient (CC), bias, root mean square error (RMSE), normalized root mean squared error (NRMSE), and Nash–Sutcliffe model efficiency coefficient (NSE). Approximately 15,000 serial model simulations (approximately two-week real-time) were needed to obtain the best parameter sets at each of the two sites using a standard, dual-core workstation. Found values of the soil parameters that did not vary spatially within the EC footprint and the ranges used during their calibration were similar to previously conducted studies [12,47,54,77] and compared favorably with results from soil pedo-transfer functions based on bulk density and particle size fractions [1,78].

During the model calibration procedure, the vegetation parameters shown in Table 6 were updated according to remotely sensed data. In most cases, the parameters utilized to characterize vegetation conditions were estimated using the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor data [12,54,86,87,88]. Albedo ($\alpha$) was obtained from the bi-hemispherical reflectance (BRDF, white sky albedo) on the visible bands [89,90] from version 6 of the MCD43A3 MODIS BRDF daily product. NDVI was obtained at a daily time scale from the MCD43A4 MODIS 500 m surface reflectance product. The other dynamic parameters described in Table 6, such as canopy field capacity (S), free throughfall coefficient (p), and optical transmission coefficient (k${}_t$), were estimated based on the 500 m MCD15A3H leaf area index (LAI) four-day composite product. The light extinction coefficient (k) used in the Lambert equation to estimate k${}_t$ at the ARM-CF site was 0.62 +/− 0.17 (croplands) and 0.5 +/− 0.15 (grasslands) for the MOISST site [80]. Stomatal resistance (r${}_s$) is based on the energy-limited relation [81,91] that depends on photosynthetic active radiation (PAR) and LAI for scaling up from the leaf to the canopy level [81,82]. r${}_s$ was estimated on an hourly time scale to reproduce the diurnal changes of vegetation, since the stomatal openings are influenced by weather conditions and changes during the day [92]. PAR is driven by SW radiation that is measured every 30 min at the EC systems and can be calculated based on the fraction of SW radiation absorbed by vegetation. Therefore, leaf-level r${}_s$ calculations were conducted through the use of SW data and the maximum seasonal stomatal conductance (g${}_{max}$) that is a biome-specific value [35]. Finally, on the assumption that vegetation cover (v${}_f$) and LAI do not experience abrupt changes during a regular day, daily LAI values were calculated using linear interpolation from the four-day LAI MCD15A3H MODIS product; daily LAI values were then used to upscale the hourly time series of leaf-level r${}_s$ to the canopy level.

A model validation was performed at each of the calibration sites but for a different time periods than the ones used during the calibration (see Table 4). During the validation stage, atmospheric forcing and dynamic vegetation parameters (that depend upon site and remote sensing data) are the only elements that changed. The validation metrics to assess the model performance were, again, the CC, bias, RMSE, NRMSE, and the NSE.

2.4. Model Transferability Approach and Evaluation

The static model parameters found through calibration (see Table 5) at ARM-CF (cropland) and MOISST (grassland) were then directly transferred to US-A74 (cropland) and US-A32 (grassland), respectively, to test the model prediction skill and parameterization robustness relative to simultaneous hourly observations of the predicted variables, during years different to the calibration or validation procedures (see Table 4). The time period for the transferability test was determined based on data availability and continuity. Thus, calibration or recalibration was not performed at this stage. This guaranteed another level of independent verification of the model performance across geographically and climatically similar sites. Although the static parameters were directly transferred, dynamic vegetation parameters at US-A74 and US-A32 were still informed by MODIS. Appraisal of the performance of the model transferability relied on the CC, bias, RMSE, NRMSE, and NSE. The results from this section provided an assessment of the model’s capacity for operating as a complementary tool for energy flux estimations (including ET), soil temperatures, and moisture at standard weather stations that measure the forcing needed to run the model.

3. Results

3.1. Model Training and Validation

The OAT precalibration procedures showed that, at both ARM-CF and MOISST, soil pore-size distribution index (m), heat capacity (C${}_s$), and saturated hydraulic conductivity (K${}_S$) were the parameters that mostly influenced the simulation results. For MOISST, besides those, the soil air-entry pressure (${\mathsf{\Psi }}_b$) and hydraulic conductivity decay (f) were also significant. Previous studies outlined the importance of these parameters during the water infiltration and shallow-soil energy exchange processes [93,94]. The parameter values found during the calibration process using the SCE algorithm are summarized in

Table 7. tRIBS static soil and vegetation (underlined) parameter values found through calibration at both ARM-CF and MOISST.

Parameter	ARM-CF	MOISST	Units
K${}_S$	21.84	4.85	[mm/hr]
${\theta }_s$	0.552	0.61	[]
${\theta }_r$	0.017	0.11	[]
m	0.57	0.52	[]
${\mathsf{\Psi }}_b$	−0.373	−99.2	[mm]
f	5.00 × 10${}^{-7}$	0.07	[mm${}^{-1}$]
A${}_s$	1.109	388	[]
A${}_u$	1.109	388	[]
n	0.431	0.51	[]
k${}_s$	0.989	1.6	[J/msK]
C${}_s$	1.061 × 10${}^6$	1.383 × 10${}^6$	[J/m${}^3$K]
K	0.2911	0.2911	[mm/hr]
b${}_2$	3.209	3.527	[mm${}^{-1}$]
H	0.2953	0.4476	[m]
${\theta }^{*}$	0.55	0.4939	[]
${\theta }_s$	0.1792	0.1577	[]

.

The time series of the remotely sensed derived v${}_f$, $\alpha$, LAI, S, p, k${}_t$, and r${}_s$ during the calibration period are illustrated in Figure 4 for ARM-CF and MOISST.

The vegetation fraction, v${}_f$, shows the largest values (i.e., V${}_f$≈ 0.7) from July to August and the lowest (i.e., v${}_f$≈ 0.1) from February to March as a result of the seasonal variations in NDVI (see Table 6) at both ARM-CF and MOISST sites. On the other hand, $\alpha$ shows both similarities and differences between sites, with values ranging from 0.20 to 0.35 year-long but without a clear pattern of seasonality. At MOISST, LAI and S are the highest between May and October as a consequence of vegetation growth. A different pattern was evident at the ARM-CF site where the maximum LAI for the calibration period occurred in April. Both p and k${}_t$ show similar temporal patterns at the same station that appear to be higher between November and May. However, ARM-CF presents higher values year-round for both variables (0.75 ≤ p ≤ 1 and 0.9 ≤ k${}_t$≤ 1) as a result of lower LAI values compared to MOISST (0.1 ≤ p ≤ 0.8 and 0.4 ≤ k${}_t$≤ 0.8), which means a less dense vegetation cover in ARM-CF year-round. Finally, r${}_s$ presents high values (r${}_s$≈ 50 s/m) from November to April at both stations, although at ARM-CF, there is a sharp decrease in r${}_s$ from mid-November to mid-December. The lowest values of r${}_s$ (r${}_s$≈ 10 s/m) occur during summertime at both stations. In general, r${}_s$, k${}_t$, and p are inversely related to S and LAI. This is the reason why the first set of parameters tend to be maximum during winter. The model results of the simulated vs. observed hourly values and statistical simulation skill metrics for the model training stage are illustrated in Figure 5, Figure 6, and

Table 8. ARM-CF and MOISST calibration metrics for each of the output variables shown in Figure 5 and Figure 6 with respect to observed values at hourly (HH) and daily (DD) resolutions. Net radiation (NR, W/m${}^2$), latent heat flux (LE, W/m${}^2$), sensible heat flux (H, W/m${}^2$), ground heat flux (G, W/m${}^2$), soil surface temperature (SST, ${}^{\circ }$C), surface soil moisture (SSM, -), and root-zone soil moisture (RSM, -). The statistical metrics are correlation coefficient (CC), bias, root mean squared error (RMSE), normalized RMSE (NRMSE), and Nash–Sutcliffe model efficiency coefficient (NSE).

Station →		ARM-CF					MOISST
Variable →		NR	LE	H	G	SST	NR	LE	H	SST	SSM	RSM
CC	HH	0.91	0.78	0.83	0.77	0.92	0.9	0.81	0.87	0.86	0.64	0.87
	DD	0.95	0.75	0.73	0.76	0.98	0.37	0.83	0.46	0.94	0.62	0.87
Bias	HH	0.24	−0.36	−8.99	−8.76	0.08	0.68	1.54	−0.35	−0.02	−0.26	−0.05
	DD	0.14	−0.36	−9.07	−8.6	−0.31	0.68	1.72	−0.34	−0.02	0.01	−0.05
RMSE	HH	77.55	66.64	60.54	22.13	2.99	126.20	63.77	85.43	5.79	0.19	0.02
	DD	24.57	39.47	43.15	9.97	2.7	78.13	42.47	31.56	3.44	0.17	0.02
NRMSE	HH	0.96	1.56	1.94	−4.11	0.17	0.98	1.12	2.27	0.30	0.47	0.10
	DD	0.31	0.92	1.38	-1.86	0.15	0.58	0.76	0.84	0.18	0.44	0.10
NSE	HH	0.81	0.56	0.59	0.36	0.84	0.55	0.66	0.71	0.83	0.20	0.68
	DD	0.85	0.43	0.27	0.23	0.94	0.22	0.17	0.48	0.96	0.18	0.77

.

From Table 8, the calibration results show that, overall, model skill is higher for hourly (HH) compared to daily averaged (DD) values. At hourly time steps, all simulated variables show NSE ≥ 0.5 (CC ≥ 0.78) except for G (NSE = 0.36, CC = 0.77) at ARM-CF and SSM (NSE = 0.2, CC = 0.64) at MOISST, with NSE > 0 in both cases. NSE values greater than 0.5 represent simulations with sufficient quality while NSE > 0, although not outstanding, provides confidence that the estimations are better than the historical time series means. Scatterplots of Figure 5 and Figure 6 also show that, generally, model simulations and observations agree well across the range of values except possibly by SSM and RSM at MOISST. Some positive biases (meaning model overestimation) are observed for LE, H, and SST at both (ARM-CF and MOISST) sites. Negative biases (i.e., model underestimations) occur for high values of NR and low values of SST also at both sites. At MOISST, simulations of SSM present large, persistent biases, most significantly at low values, although those persistent errors are also present across the entire distribution of hourly values. RSM presents a similar behavior to SSM at MOISST, although the biases appear less pronounced. Although the simulation skills at the daily (DD) time step illustrate higher scores for NR and SST at ARM-CF and for SST and RSM at MOISST, they also show lower scores, particularly for LE, H, and G at ARM-CF and NR, LE, H, and SSM at MOISST. Nonetheless, NSE for daily (DD) comparisons shows that the simulations are still better than the historical means in all cases.

The daily-aggregated time series of observed precipitation and predicted variables at ARM-CF and MOISST for the best calibration parameter sets are displayed on Figure 7 and Figure 8, respectively, during the entire calibration period. The time series are accompanied by envelopes illustrating the daily ranges of variability through their standard deviation. Due to its high hourly persistence, SSM and RSM’s variational envelopes are particularly narrow. At both calibration sites, the seasonal patterns of fractional distribution of the incoming solar radiation into LE, H, and G (and thus NR) are well captured by the model simulations and their daily variability closely follows. Near zero and negative NR (meaning net energy deficits) tend to occur most frequently during the winter season but positive values dominate the warm months. LE tends to be positive year-round except during the short times of water vapor condensation. LE is highest during the summer and spring seasons. H presents average positive values (meaning upward heat transfer) during the warm months but zero to negative (i.e., downward transfer) during winter. G has fluctuating (i.e., positive and negative) values year-round, meaning quickly changing patterns of downward and upward soil heat transfer in short periods of time. Overall, at both stations, LE reaches the largest peak values (∼400 W/m${}^2$) while G presents the smallest absolute values (∼50 W/m${}^2$) among all heat fluxes. SST closely follows the seasonal changes in incoming solar radiation with peak average values around 40 ${}^{\circ }$C in summer and near (i.e., slightly above or below) 0 ${}^{\circ }$C during winter. SSM and RSM do not appear to show a marked seasonal pattern. Instead, P plays a significant role in increasing both SSM and RSM and simultaneously reducing SST and H. P also seems to increase LE during the warm season. SSM is the poorest simulated variable at MOISST, particularly when P is high. Across sites and overall (except possibly by SSM at MOISST) predicted seasonal variability including average, maximum, and minimum values indicate model accuracy and appropriateness.

Figure 9 and Figure 10 provide a closer look at the sub-daily cycles of the variables of interest for a ten-day summer period but during different years (2004 for ARM-CF and 2014 for MOISST) including pre-and post-storm periods. Both figures indicate that simulations capture the surface energy partitioning and soil temperature especially for ARM-CF. The presence of significant rainfall events (particularly those occurred in ARM-CF) flattens such diurnal cycles for a short time, but then, the cycle pattern repeats as in pre-storm conditions. During this particular time, simulations at MOISST lacked accuracy for NR, LE, and of course the persistent positive biases in SSM and RSM. Analogously to Figure 8, SSM and RSM do not seem to follow a diurnal cycle but are heavily driven by precipitation events. tRIBS appears to represent those sub-daily changes in shallow and root-zone soil moisture by capturing peaks and the raising and falling limbs.

The model validation results for ARM-CF and MOISST are summarized in

Table 9. ARM-CF and MOISST validation metrics with respect to observed values at hourly (HH) and daily (DD) resolutions. Net radiation (NR, W/m${}^2$), latent heat flux (LE, W/m${}^2$), sensible heat flux (H, W/m${}^2$), ground heat flux (G, W/m${}^2$), soil surface temperature (SST, ${}^{\circ }$C), surface soil moisture (SSM, -), and root-zone soil moisture (RSM, -). The statistical metrics are correlation coefficient (CC), bias, root mean squared error (RMSE), normalized RMSE (NRMSE), and Nash–Sutcliffe model efficiency coefficient (NSE).

Station →		ARM-CF					MOISST
Variable →		NR	LE	H	G	SST	NR	LE	H	SST	SSM	RSM
CC	HH	0.92	0.74	0.78	0.74	0.88	0.9	0.82	0.82	0.85	0.73	0.62
	DD	0.67	0.55	0.59	0.54	0.91	0.60	0.79	0.19	0.95	0.74	0.62
Bias	HH	0.13	−0.32	−16.55	−26.45	0	1	0.41	0.17	−0.06	−0.15	−0.17
	DD	−0.10	−0.30	2.61	3.04	−0.02	1	0.52	0.17	−0.03	−0.14	0.24
RMSE	HH	78.47	63.39	71.97	20.86	6.30	130.5	46.25	69.45	5.45	0.16	0.093
	DD	62.35	41.75	45.32	9.98	3.02	78.35	28.17	30.01	2.85	0.15	0.079
NRMSE	HH	0.94	1.45	2.15	−2.40	0.41	1.13	0.95	2.17	0.29	0.28	0.21
	DD	0.61	0.87	1.09	−1.17	0.13	0.68	0.59	0.96	0.15	0.27	0.21
NSE	HH	0.83	0.51	0.53	0.50	0.74	0.45	0.34	0.71	0.81	0.55	−0.10
	DD	0.02	0.23	0.20	0.09	0.81	0.08	0.32	−0.17	0.93	0.55	−0.10

. Consistent with the calibration assessment results, model simulation skill for hourly (HH) was higher than for daily (DD) data during the validation period. At hourly time resolution, such results illustrate NSE ≥ 0.5 for all predicted variables at ARM-CF (including G). However, MOISST showed lower values overall except for SSM and SST than ARM-CF yet greater-than-zero NSE scores. An exceptionally low value (NSE = −0.1) was obtained for RSM at MOISST that contrasts with the high value (NSE = 0.68) during calibration. Similar to the calibration results, the seasonal and sub-daily time series of simulated variables mimic the observed values.

Taking calibration and validation together, at ARM-CF, the variables that are better simulated by the model are, in order, NR (NSE = {0.81, 0.83}, CC = {0.91, 0.92}), SST (NSE = {0.84, 0.74}, CC = {0.92, 0.88}), H (NSE = {0.59, 0.53}, CC = {0.83, 0.78}), LE (NSE = {0.56, 0.51}, CC = {0.78, 0.74}), and G (NSE = {0.36, 0.50}, CC = {0.77, 0.74}). At MOISST, the best-simulation rank is SST (NSE = {0.83, 0.81}, CC = {0.86, 0.85}), H (NSE = {0.71, 0.71}, CC = {0.87, 0.80}), NR (NSE = {0.55, 0.45}, CC = {0.90, 0.90}), LE (NSE = {0.66, 0.34}, CC = {0.81, 0.80}), and SSM (NSE = {0.20, 0.55}, CC = {0.64, 0.73}).

3.2. Model Inter-Site Transferability

The results of the model parameter transfer experiments are shown in Figure 11, Figure 12, and

Table 10. ARM-CF to ARM-A74 and MOISST to ARM-A32 model parameter transfer evaluation statistics with respect to observed values at hourly (HH) and daily (DD) resolutions. Net radiation (NR, W/m${}^2$), latent heat flux (LE, W/m${}^2$), sensible heat flux (H, W/m${}^2$), ground heat flux (G, W/m${}^2$), soil surface temperature (SST, ${}^{\circ }$C), root-zone soil temperature (RST, ${}^{\circ }$C), and surface soil moisture (SSM, -). The statistical metrics are correlation coefficient (CC), bias, root mean squared error (RMSE), normalized RMSE (NRMSE), and Nash–Sutcliffe model efficiency coefficient (NSE).

Station →		ARM-A74						ARM-A32
Variable →		NR	LE	H	G	SST	RST	NR	LE	H	G	SST	RST	SSM
CC	HH	0.92	0.69	0.70	0.76	0.79	0.98	0.88	0.69	0.73	0.59	0.80	0.99	0.59
	DD	0.70	0.51	0.06	0.83	0.96	0.98	0.58	0.49	0.16	0.77	0.96	0.99	0.58
Bias	HH	0.13	0.13	−0.32	13.27	−0.05	−0.07	−0.12	−0.15	−0.68	−4.10	−0.10	−0.10	0.80
	DD	0.03	0.71	0.11	13.56	−0.05	−0.06	−0.11	0.08	−0.59	−0.07	−0.14	−0.12	0.71
RMSE	HH	78.47	92.60	106.90	33.20	6.70	1.99	105.61	75.99	122.02	28.80	8.852	3.81	0.16
	DD	38.59	76.59	65.79	9.61	2.46	1.96	59.66	50.03	74.26	5.66	4.61	4.02	0.15
NRMSE	HH	0.76	−.93	2.30	5.33	0.34	0.10	1.04	0.95	4.11	25.36	0.49	0.22	0.53
	DD	0.34	0.75	1.26	1.54	0.12	0.10	0.54	0.57	2.10	3.83	0.25	0.22	0.49
NSE	HH	0.83	0.31	0.56	−0.11	0.61	0.93	0.77	0.60	0.60	0.36	0.60	0.87	0.18
	DD	0.37	−0.68	−1.33	0.08	0.91	0.93	0.31	0.45	−0.01	0.57	0.81	0.86	0.16

for ARM-A74 and ARM-A32, respectively.

At the hourly time resolution, the set of statistical values from Table 10 ensure high (NSE ≥ 0.5) to appropriate (NSE ≥ 0) simulation skill in most of the target variables, in consistency with the results at ARM-CF and MOISST. For example, at ARM-A74 and ARM-A32 (Table 10), although not calibrated, RST showed exceptionally high scores (NSE = {0.93, 0.87}, CC = {0.98, 0.99}). Similarly, NR (NSE = {0.83, 0.77}, CC = {0.92, 0.88}), SST (NSE = {0.61, 0.60}, CC = {0.79, 0.80}), and H (NSE = {0.56, 0.60}, CC = {0.70, 0.73}) had all NSE ≥ 0.5, meaning great skill in capturing the temporal dynamics of those variables. The predictions of LE, however, decreased in quality at ARM-A74 (NSE = 0.31, CC = 0.69) compared with the ones obtained during the calibration and validation procedures (NSE ≥ 0.5) at ARM-CF. Nonetheless, at ARM-A32, LE simulations (NSE = 0.60, CC = 0.69) maintained consistent quality with the results during the calibration process at MOISST (NSE = 0.66, CC = 0.81). Lastly, simulations of G showed a poor performance (NSE = −0.11, CC = 0.76) at ARM-A74 and only acceptable scores at ARM-A32 (NSE = 0.36, CC = 0.59). Consistent with the model assessment during calibration and validation, the skill at daily time steps showed improvements for SST, RST, and G but decreased performance in all other variables. In some cases, the transferability of the model for evaluation at the aggregated daily time steps is not even recommended for use. This is the the case for LE and H at ARM-A74 and for H at ARM-A32.

4. Discussion

The main objectives of this study were (1) to evaluate the capability of a physical based land-surface model that assimilates remote sensing data on vegetation dynamics to simulate surface energy fluxes, surface and root-zone soil temperature, and moisture at the scale of the eddy flux footprint and (2) to test the seamless transferability of this tool at sites with similar soil and vegetation characteristics and standard weather station measurements. The overall goal was the development and testing of a modeling platform that could potentially be used as real-time virtual eddy covariance tower and soil monitoring station in ungauged sites where this information is not measured (e.g., standard weather stations). One of the benefits of using a process-based model is its explicit declaration of internal states, robustness, and applicability to regions where accurate input information (e.g., soils, vegetation, and atmospheric forcing) is available. However, this is also precisely one of the limitations too, as this information is not always available. The multi-physics framework is supported through equations of the processes occurring at each site to obtain the right answers for the right reasons. One of the key aspects of correctly representing the energy- and water-budget-related processes is the consideration of vegetation dynamics (i.e., canopy cover and activity). Contrary to the great majority of previously published remote-sensing based models [25,26,27,28,29,30,31,32], tRIBS successfully assimilates LAI, NDVI, and albedo from satellite sources, and ancillary soil and terrain data so that land surface-related parameters dynamically evolve according to physics-based equations that determine important canopy and ground hydraulic and thermal processes such as solar energy reflection and absorption, incoming and outgoing short- and longwave radiation, vegetation shading, atmospheric and soil thermal conductivity, soil heat diffusion and temperature damping depth, precipitation water interception and its evaporation, bare soil evaporation, plant transpiration, shallow water infiltration and motion of wetting fronts, runoff generation, and subsurface water flow, among others.

Model training and validation proved the robustness of the found static and dynamically changing parameters and provided confidence in the model performance at the same sites where calibration was conducted. Nonetheless, since the model was calibrated with hourly data, its performance was, overall, better than when it was assessed against daily-aggregated time series. This could indicate that the inclusion of daily-aggregated targets in the calibration procedure could ameliorate the cumulative biases introduced by only using hourly data. One of the key learned lessons from model training and validation is that the level of interconnectedness of processes requires high-quality training datasets. Otherwise, errors in at least one training set will be propagated into errors (that could be multiplicative) in other dependent variables. For instance, soil moisture is a key variable driving the magnitude of almost all other predictors in this study. Therefore, obtaining high-quality field measurements of soil water content is a sine qua non condition for conducting an unbiased model training. In our case, it was difficult to find high-quality observations of soil moisture. We could attribute much of our uncertainties to inaccurate to poor-quality soil moisture data.

The model inter-site transferability experiments provided another degree of model validation at not only a different time but also a different (paired) location with some ecosystem similitude. The ability to perform a straight-forward model transfer, relying on vegetation dynamics, becomes a paradigm shift in the simulation framework of land surface models that overrelied on exhaustive model calibration procedures of fixed vegetation parameters for each work-site. Overall, the simulated hourly time series demonstrated the existence of parameter sets that can be transferred to similar ecosystem locations to rapidly obtained ballpark estimations in ungauged (i.e., without eddy covariance measurement) sites. The model estimates were better than the historical means while avoiding equifinality issues that usually hinder model parameter transfer in space and/or time. Future research should focus on improvements of this methodology, perhaps facilitating the modeling of multiple locations at once.

Among the limitations of this study, we can mention (1) the uncertainties of the weather forcing measurements, (2) errors in the observed time series of the modeled variables, (3) the space and time resolution of the satellite imagery, and (4) the low number of test sites. First, inaccurate weather forcing negatively influences the performance of any land-surface or hydrologic model as predictions heavily depend on first-order variables such as incident solar radiation (SW) or precipitation (P) [95]. For instance, after a careful inspection, we found that, at ARM-CF, SW had a significant number of gaps across all years. Therefore, we deferred to FLUXNET data that are gap-filled using a synthetic method [22]. At MOISST, despite measurements of SW being inspected, tower managers perform their own SW gap-correction method without subsequent quality-control assessment. This could be one of the reasons why NR does not reach high simulation skill scores at this site. Second, there were some constraints in finding reliable, continuous, and depth-matching soil moisture, soil temperature, and ground heat flux data. For instance, SSM time series were not available or showed poor quality when available at ARM-CF and ARM-A74 and, therefore, we decided not to trust them. Likewise, RSM was only available at MOISST, given that this is a long-term soil moisture intercomparison benchmark site. RST was only measured at ARM-A74 and ARM-A32 but not at the sites used for calibration and validation. Other issues that possibly impacted the quality of the SSM and RSM simulations are the significant discrepancies between soil moisture sensors of different brands [50,96] in addition to their poor performance during below-zero temperatures [97,98] including the ground heat plates, which produce an error in the ratio of downward to upward flux of shortwave radiant energy estimations [99]. Additionally, matching the depth at which SSM, RSM, SST, and RST were measured and simulated becomes an issue as tRIBS outputs these values at 10 cm and 100 cm depths. The majority of the eddy towers had soil observations at 2.5, 5, 60, and 90 cm below ground. Although we tried to match the closest depths, we believe this might have slightly contributed to the modeled and observed discrepancies. The errors in the estimation of LE can be explained, in part, by the mixed quality of estimation of the soil water content (i.e., SSM and RSM). Finally, the time series of G, as the lowest energy flux, showed numerous gaps and quality flags that hindered a fair comparison with our model simulations. The quality of the flux heat plates had previously been questioned and discussed in [100]. Third, the limiting resolution of the dynamic vegetation-related information (from MODIS) hindered the model capacity to spatially and temporarily cope with sub-pixel and sub-daily phenologic changes that might occur, for instance, after summer precipitation events. Considering 500 m aggregated, instead of data at the resolution of the tower foot-print (i.e., smaller than 200 m) could result in sensing errors of the vegetation dynamics or in accounting for processes that are not necessarily occurring within the foot-print of interest. Future work might be needed to assess the model performance when using higher-resolution products such as Landsat 30 m or Sentinel 20 m to estimate the vegetation components. Another well-known limitation of satellite images is the presence of clouds that obscure ground status. This issue, which is particularly critical during winter and spring in Oklahoma, was overcome through image interpolation between two or three successive images. Finally, the lack of hourly temporal resolution from the satellite images inhibits replicating diurnal cycles of vegetation activity. One of the critical vegetation parameters that control the rate of transpiration is the stomatal resistance (r${}_s$). This parameter was estimated at an hourly time step by using LAI and FPAR four-day composites that were linearly interpolated to hourly by using time series of SW. This estimation could have been more precise through direct use of the photosynthetic photon flux density (PPFD) data from the sensor typically mounted in many EC towers. This way, a direct estimation of FPAR could have provided a more accurate estimation of the plant-use of solar energy at the foot-print scale. However, since this study aimed to test the model at places without an eddy covariance system, obtaining this parameter through satellite imagery would be the best (but not available yet) alternative. Fourth, while the results provide confidence in the potential ability of the modeling framework to successfully perform across similar ecosystems, the tests conducted in this study were limited to grassland and cropland environments within the U.S Great Plains. Without doubt, more testing is needed including different ecosystems and climates.

Future work might include the development of a stakeholder-oriented virtual tool where tRIBS can be run online and at any point of interest with weather forcing data (e.g., Mesonet or re-analysis model outputs). Moreover, it could be utilized in conjunction with drone or sub-orbital technologies that would refine the remote sensing images to the tower footprint, which would help improve the spatial resolution of vegetation parameters (such as r${}_s$) and, thus, the overall simulation accuracy of the results.

5. Summary and Conclusions

This study compiled atmospheric, soils, vegetation, and hydrologic information from four eddy-covariance systems within the U.S. Great Plains to evaluate the capabilities of a process-based model that uses standard weather station measurements and remote sensing of vegetation to estimate the components of the surface energy budget, soil temperatures, and water content at the surface and the root zones. The results of the model training, validation and parameter transferability tests allowed us to make the following conclusions:

1.

Hourly simulations adequately predicted NR, SST, RST, and H, including the representation of seasonal variability and daily cycles. Even though LE and SSM sometimes, obtained mixed skill assessments, the simulations always guaranteed a quality level better than historical means and, in some cases, with sufficient quality.

2.

Model validation proved the robustness of the found static and dynamically changing parameters and provided confidence for the model performance at the same sites where calibration was conducted.

3.

Inter-site transference of the model framework showed that the model was able to assimilate data about vegetation dynamics from the remotely sensed information to update important in situ vegetation parameters that resulted in adequate hourly model performance metrics. Model transferability experiments from ARM-CF to ARM-A74 (cropland) and from MOISST to ARM-A32 (grassland) provided arguments to explore future use of a set of precalibrated static parameters with another set of dynamically evolving, in situ, vegetation parameters between regions of pedologic, ecosystem, and hydrologic similarity.

4.

The straightforward model transfer relying on vegetation dynamics marks a paradigm shift in the simulation framework of land surface models that has overrelied on exhausting model calibration procedures of fixed vegetation parameters for each work-site. Parameter sets that can be transferred to similar ecosystem locations become powerful tools for facilitating the modeling of multiple locations at once.