Evapotranspiration in Northern Agro-Ecosystems: Numerical Simulation and Experimental Comparison

Evapotranspiration and near-surface soil moisture dynamics are key-entangled variables regulating flux at the surface-atmosphere interface. Both are central in improving mass and energy balances in agro ecosystems. However, under the extreme conditions of high-latitude soils and weather pattern variability, the implementation of such cou pled liquid and vapor phase numerical simulation remain to be tested. We consider the nonisothermal solution of the vapor flux equation that accounts for the thermally driven water vapor transport and phase changes. Fully coupled flux model outputs are compared and contrasted against field measurements of soil temperature, heat flux, water content, and evaporation in a subarctic agroecosystem in Alaska . Two well- defined hydro-meteorological situations were selected: dry and wet periods. Numerical simulation was forced by time series of incoming global solar radiation and atmo spheric surface layer thermodynamic parameters: surface wind speed, ambient tem - perature, relative humidity, precipitation, and soil temperature and soil moisture. In this simulation, soil parameters changing in depth and time are considered as dynami - cally adjusted boundary conditions for solving the set of coupled differential equations. Results from this evaluation give good correlation of modeled and observed data in net radiation ( R net ) ( R 2 of 0.92, root mean square error (RMSE) of 45 W m −2 ), latent heat (0.70, RMSE of 53 W m −2 ), and sensible heat ( R 2 = 0.63, RMSE = 32 W m −2 ) during the dry period. On the other hand, a poor agreement was obtained in the radiative fluxes and turbulent fluxes during the wet period due to the lack of representation in the radiation field and differences in soil dynamics across the landscape.


Introduction
Northern latitudes have been identified as a region where global climate change will have earlier and stronger impacts than in other regions of the world [1][2][3][4]. Most of the region is underlain by discontinuous permafrost or perennially frozen ground in which temperatures remain below 0°C for at least two consecutive years. An active layer on top of the permafrost experiences seasonal thaws and is the primary dominant subsurface component of the land-atmosphere system [5]. Under climate warming scenario, much of this terrain would be vulnerable to subsidence, particularly in ice-rich areas of relatively warm, discontinuous permafrost, and shrinking ponds and lakes [3,[6][7][8][9][10]. All these changes will potentially alter the exchange of surface energy, water, and carbon cycles in high-latitude ecosystems [11] and consequently, the response at regional level to the atmosphere system. Soil moisture plays a critical role in the surface energy balance and water cycle in these regions [1,12,13]. It is widely recognized that the soil moisture confined in a thin layer underneath the land surface influences the partitioning of the surface energy fluxes simultaneously modifying surface thermal conductance and rates of evaporation [14]. An example of such an important role is the control of precipitation transfer into the soil and the partitioning of incoming solar radiation into latent, sensible, and ground heat fluxes [15,16]. In addition, soil moisture and temperature status affect biological processes such as soil microbial activity, seed germination, and plant growth. These variables in turn also affect water and nutrition absorption and solute transport in soil. In a climate change scenario, high-latitudes soils will experience increased summer dryness as climate warming progresses, changing therefore atmospheric vapor pressure conditions and thereby enhancing evapotranspiration (ET) rate. In terms of seasonal effects, inadequate snowmelt infiltration or rainfall during spring and early summer often causes crop water stress and reduction in yield of small grains [17,18] in agricultural activities of subarctic regions. Therefore, understanding the variation of evapotranspiration (ET) and its impact on crop growth becomes of absolute importance because it mainly controls the available soil water and, therefore, is a limiting factor in agriculture productivity and sustainability. As a result, continuous monitoring of soil water content and soil temperature is a priority in the fields of agronomy and hydrology [19].
Several modeling studies have focused on soil carbon reservoirs (e.g. [20][21][22]) and permafrost degradation in natural ecosystems across the circumpolar region [21,23]; nevertheless, agroecosystem has not been taken into consideration until a recent study by Ruairuen et al. [24]. Despite the mentioned complexities in the soil medium, similarities between high latitude and mid-latitude agricultural soils exist mainly during the growing season. This allows for making use of models that are currently in use for mid-latitude agricultural settings. In this case, a fully coupled differential equation system considering both soil temperature and vertical soil moisture distribution, and their interactions are utilized to bring emphasis on the sub-medium transport in contrast to most large-scale ecosystem models where one or two soil layers are used to simulate soil moisture dynamics in ecosystem models (e.g. [25]).
In this study, we use the numerical model developed by Bittelli et al. [26], which fully couples heat and water transport to deduce the coupled water and heat transport across the soil medium forced by radiation and meteorological conditions. As demonstrated in mid-latitudes studies by Bittelli et al. [26], this approach enables numerically stable, energy and mass conservation equation solution in terms of the external forcing and boundary conditions. Such an approach requires a modeling framework that incorporates the interactions among meteorological variables (e.g., air temperature, relative humidity, precipitation, and solar radiation) and soil properties (e.g., soil temperature, soil moisture, and soil water potential) into the coupled numerical model. We also conducted a field experiment to measure net radiation, sensible heat flux, soil moisture, and soil temperature, and use those measured parameters to calculate evapotranspiration in order to compare with simulated results.

Field experiment
The field study was conducted at the Fairbanks Experiment Farm (FEF) of the University of Alaska Fairbanks (UAF) Agricultural and Forestry Experiment Station (AFES) Fairbanks Alaska (Figure 1), USA (64°51′16.6″ N, 147°51′36.4″ W, 150 m above sea level) during summer 2013. The soil within the lysimeter plots, established in a previous study, was used for this study because large amounts of data (i.e., soil moisture, soil temperature, and soil moisture potential) were available. The soil composition was sandy loam with 66 sand, 29 silt, and 5% clay, and with the available water holding capacity of about 0.18-0.36 m 3 m −3 that was determined from a soil moisture characteristic curve [24]. Parameters for soil hydraulic properties to be introduced in the numerical simulation are listed in Table 1. Volumetric soil moisture content was measured in situ using three soil moisture sensor (10HS; Decagon Devices Inc., Pullman, WA, USA) at 5, 10, and 20 cm. The sensor has 14.5 cm long prongs, 3.3 cm wide, so the 5-cm sensor spanned the depth 3.3-6.7 cm, the 10-cm sensor spanned the depth 8.3-11.3 cm, and the 20-cm sensor spanned the depth 18.3-21.3 cm. Gravimetric samples were also taken to calibrate the 10HS sensor [24]. Soil moisture was continuously measured with an interval of 30 min. The bulk density was measured in the tested site ( Table 1). The soil temperature (S-TMB-M006, Onset HOBO Corporation, Bourne, MA, USA) was measured at depths of 5, 10, and 20 cm below the soil surface and used to obtain the ground heat flux in this study.
The observation-based meteorological parameters included air temperature (Tair), relative humidity (RH), air pressure, wind speed (u), and direction, and precipitation at 2 m height above the ground were obtained at 1-min intervals at the experimental station. One-min recordings of these data were averaged to obtain hourly data for input to the simulation.
An independent measure of evapotranspiration was determined using Penman-Monteith method and the more continuous series of data available on this period. Sensible heat flux was measured locally (i.e., ecosystem scale) by means of and eddy-covariance (EC) instrument and processed considering signal distortions under all weather conditions [27] and, at landscape scale, based on a large aperture scintillometer (LAS) [28][29][30].

Model description (PSP_coupled)
The numerical model was coded in Python and is set in a time-evolving one-dimensional simulation of coupled flow of liquid water, heat, and water vapor. The model description can be found in Bittelli et al. [26,31]. Figure 2 shows the models' conceptual scheme indicating the coupled layers and driving boundary conditions, i.e., soil temperature, liquid water

Soil property Value
Bulk density (g cm −3 ) 0.7 Air entry potential (J kg −1 ) −1.5 Mass sand (kg kg −1 ) 0.66 Mass clay (kg kg −1 ) 0.05 The main file is main.py, which contains the calls to other embedded subroutines listed. The module PSP_boundary defines the initial and boundary conditions. The PSP_public contains all variables that are read by all modules such as latitude, longitude, altitude, albedo, atmospheric pressure and clay content, initial soil temperature, and soil matric potential. The PSP_soil is written to define the soil properties. The PSP_couple1D is the module that implements the solver for the different flux equations. The PSP_longWaveRadiation is for computing the long wave radiation component of the radiation balance at the soil surface. The PSP_grid module is for building the computational grid and PSP_ThomasAlgorithm for solving the system of equations. The PSP_plot and PSP_plotEnergy are modules for visualizing the data input and output from the model.

Initial setting for model simulation
The initial conditions for dry and wet periods were calculated using in situ data. To implement the simulation scenarios, two data files need to be created "soil.txt" and "weather.dat." The soil.txt file is required for data input of soil properties such as soil depth (the model set up from 0 to 1.5 m), the saturated soil moisture, residual water content, soil hydraulic properties, and soil matric potential ( Table 1). The soil file can be used to simulate the two periods; however, additional settings indicating initial conditions should be modified in the PSP_public module. The weather file is required at hourly weather parameters such as solar radiation, air Figure 2. Scheme of the computational grid with the driving force terms (temperature, soil water potential, and soil vapor concentration), the soil conductivities and the resistances involved at the soil-atmosphere interface (adapted from Bittelli et al. [26]). The black dot is the mass balance for heat flow and water flux at a given node.
temperature, precipitation, relative humidity, and wind speed as an input to calculate related parameters. Time step is set to 300 s and input data of 1-h resolution.
The mass of sand, silt, clay, and bulk density was obtained from in situ measurements. The parameter for the hydraulic properties was obtained from Cambell and Shiozawa [32], K s is saturated hydraulic conductivity, θ s is the saturated soil moisture content, and b is constant value.
The PSP_public is the file that needs to be adapted in all parameters that are read by all modules for the given area in which the simulation is carried out. In this case, the FEF site-specific information was input including latitude, longitude, and altitude. Moreover, the soil initial conditions such as soil water potential, soil temperature, and albedo for dry and wet scenarios needed to be applied into this module ( Table 2).
In this study, the value of albedo was set as 0.2 for the dry period [33], while a value of 0.15 was applied for the wet period in agriculture land in subarctic region according to previous studies in the same agricultural setting [34].

Simulation scenarios
As mentioned previously, the model was applied to two selected periods (dry and wet) in an agricultural land described in Section 2 during the summer 2013, and its performance was evaluated. In addition and based on local meteorological information, it was verified that the atmospheric boundary layer developed forced by surface and near surface flow conditions without presence of multilayered thermal inversions [35]. However, this condition is difficult to maintain when precipitation arises. The dry period (no precipitation event) spun from 26 July to August 3 (Julian day 207-216) and wet period from 25 to 30 August (Julian day 237-242).

Dry period
The experimental data were taken from the plots that monitored soil moisture and soil temperature. The meteorological parameters measured about 10 m away from the plot for the dry period are given in Table 3 and Figure 3a.

Wet period
The meteorological conditions during the wet period 18-30 August 2013 (Julian 230-242) was cooler than the dry period in terms of an average hourly air temperature and soil temperature (Figure 3b), while there was a slight difference in solar radiation compared to the dry period (see Table 3). The RH was approximately 77% with low level of VPD (0.36 kPa) on average during the wet period ( Table 3). A total precipitation of 37.60 mm was also reported in this period. However, only data during 25-30 August 2013 (Julian 237-242) are used as the wet period for the simulation in this study.

Net radiation
Net radiation (R net ) is the main input in the energy balance equation driving evapotranspiration process. Giving this importance, a comparison between the net radiation observed (R net obs) and modeled (R net mod) was performed. Figure 4a shows simulated and measured hourly data of net radiation as a function of time for the dry period (26 July-3 August 2013). Results show a good agreement between modeled and observed R net with the correlation coefficient (R 2 ) of 0.92 and root mean square error (RMSE) of 45 W m −2 (Figure 5a) during the dry period. There are about 6 days out of nine total days that the model performed remarkably well to simulate the R net and the maximum difference did not exceeded 157 W m −2 on 27 July (Julian day 208).
Overall, R net mod overestimated R net obs by about 16%. Similar results have been reported by Ortega-Farias et al. [36] in a commercial vineyard where the R net simulation shows a good agreement with the measured R net (R 2 = 0.92) with the RMSE less than 48 W m −2 . On the other hand, a lower correlation between R net obs and R net mod was found during the wet period with the R 2 of 0.72 and RMSE of 67 W m −2 (Figure 5b). However, there were two days from 25 to 26 August that the model simulated very well compared with the measured R net , while it was underestimated in the following 2 days (27-28 August) and overestimated in the last 2 days during the wet period (Figure 4b). The maximum underestimation for the entire time series was about 85 W m −2 , and the maximum overestimate was 199 W m −2 .
The high and low correlation between measured and modeled R net in different periods can be related to the cloudiness variability in high latitudes. This may be explained based on the variability of incoming solar radiation between those periods as we can see in Figure 6. An average of solar radiation of 215 W m −2 with the maximum and minimum of 680, 150 W m −2 respectively were reported in the dry period, while a lower mean and maximum R net were found in the wet period ( Table 3). Ortega-Farias et al. [37] indicated that errors in the calculation of R net over a well-irrigated festuca grass were associated with the estimation of atmospheric radiation under cloudy sky conditions. Therefore, to summarize this point, R net mod was able to estimate net radiation with a good degree of precision during the dry period than during the wet period.    from the vsimulation was 15 mm with the mean of 1.67 mm d −1 for the dry period. In contrast, during the wet period, the cumulative ET from the observed was about 6 mm over 6 days, whereas only about 2 mm was found from the simulation of ET cycles. An average precipitation observed was 1.0 mm d −1 , while only 0.3 mm d −1 was found for the simulation during the wet period.

Ground heat fluxes
The average daytime ground heat flux (G) contribution to vapor flux was in the order of 26 W m −2 (ranging from 1 to 79 W m −2 ), while the value output from the simulation was on average 25 W m −2 (ranging from 0.8 to 49 W m −2 ; Figure 8). In Figure 8a, modeled G resulted as overestimated for the first two days, and then closely approached the G observed at day 3, and underestimated the observations during the last 2 days of the dry period. The maximum reported difference was 44 W m −2 during the overestimating period, whereas the same values were 32 W m −2 found during the underestimating period. On the contrary, the simulation of G during the wet period was overestimated for the entire period (Figure 8b) with the maximum difference of more than 100 W m −2 during daytime.

Sensible heat
The sensible heat (H) fluxes obtained from measurements and from the simulation for agricultural field were compared during the dry and wet period. The time series composite of hourly diurnal variation during dry and wet was illustrated in Figure 9. In the dry period, results showed a good correlation between H measured by EC and simulated from models with R 2 = 0.63 and RMSE = 32 W m −2 , whereas a lower correlation was found in H observed by LAS (R 2 of 0.52 and RMSE of 40 W m −2 ) (Figure 9a). Nevertheless, LAS measurements represent a larger fluxing footprint area than EC-derived fluxes, and the simulations output converge well between these two scale-dependent observations. Overall simulated outputs overestimated H from EC except

Soil temperature and soil moisture
Soil temperature was measured at the same depth as soil moisture. Figure 10 shows simulated and observed soil temperature as a function of time during dry and wet periods. During the dry period, the value of soil temperature from experiment reached 17°C and was higher than the simulated 11°C for the same depth at the beginning of the simulation period. However, after midday in the first day, the simulated value was higher than observed soil temperature over entire dry period (Figure 10a). The soil temperature was better predicted by the numerical model during the wet period, with an R 2 of 0.59 and RMSE of 1.82°C than the dry period (R 2 = 0.47 and RMSE = 3.27°C). The maximum difference between observed and simulated soil temperature was about 8°C during mid of the day in dry period and the times series followed each other. There was only a day that the model underestimated in a wet period, while most of the time, simulated values overestimated the observed data (Figure 10b). A large difference of the soil temperature can cause differences in the land-atmosphere temperature gradient that affect the ground heat flux as described in the previous sections. However, in this case, we found the model in good agreement in predicting soil temperatures for this environment.
The soil moisture was measured in the plot at three depths. From the measurements, we found that the high moisture content was in the lower depth than in the surface layer. The initial soil moisture during the dry period above the soil surface was about 0.28 m 3 m −3 (5 cm depth). The numerical model gave a lower level soil moisture around the same depth with a large difference of 0.21 m 3 m −3 when compared to the observed data. Some disagreement between modeled and observed data was also found during the wet period where the measurement of soil moisture from the plot was 0.30 m 3 m −3 for the first day; however, the simulation gave a lower value of soil moisture with a difference larger than 0.13 m 3 m −3 . In addition, an underestimation of simulated soil moisture potential is also reported in this study during both periods. This could be due to the lack of representation on the hydraulic properties of the soils especially for the subarctic soil. This is important because the hydraulic conductivity versus the soil moisture potential curve is highly nonlinear and, therefore, the flow of soil moisture from the upper layer to the lower layer in a wet period leads to a large decrease in hydraulic conductivity and liquid water distribution [26,31], while during the dry period, the soil moisture was more constant along the depths with less dependence on the liquid fraction. The soil moisture content fluctuated through the day according to the vapor flux as reported previously [38][39][40]. As such, improvements in the model's representation of both soil moisture and soil moisture potential in order to have an optimal simulation output. There is also a need to further study the vertical soil properties along the soil depth under agricultural land in subarctic region. Evaluation of simulation performance for subarctic soil and weather seems that more parameters might be needed to improve model simulation because of influence of the permafrost, soil properties across landscape and weather variability.

Conclusions
An effort was undertaken to simulate fluxes from surface-atmosphere exchanges based on numerically solving the coupled vapor and liquid water differential equations prescribing soil properties and turbulent exchange parameters. Numerical simulation was forced by meteorological data, radiation, and precipitation from a high-latitude agricultural farm. Similarly,  dynamic boundary conditions were introduced throughout the simulation including soil temperature, soil moisture, and soil hydraulic properties.
After examining simulation outputs and comparing them to collocated micrometeorological data, it can be concluded that time series of fluxes during the dry period seemed to be reproduced fairly better than those obtained during the wet period. In general, R net has a good agreement between modeled and observed data for both periods with RMSE of 45 W m −2 in a dry period and 48 W m −2 during a wet period. The LE also was well predicted by the model with RMSE not exceeding 53 W m −2 . This difference in turbulent fluxes agrees well with other studies in the same area over highly heterogeneous terrain with further canopy complexity [41,42]. On the other hand, G was overestimated with the maximum difference of more than 100 W m −2 in a wet period and 44 W m −2 in a dry period. In addition, the measured H by EC and LAS instruments correlated well with the model in terms of the RMSE being in the range of 32-40 W m −2 which falls within the interval of fluxing difference across landscapes observed on the same area for heterogeneous surfaces [41,42]. However, only the soil temperature correlated better during a wet period than a dry period, while the soil moisture and soil moisture potential was underestimated compared to the observed values. The low correlation in the wet period was due to significant influences of the synoptic variability introducing large changes in cloud coverage and precipitation that are difficult to reproduce by a single point one-dimensional model formulation. Nevertheless, dry conditions that are by far the most stringent conditions for agriculture sustainability reproduces well.
There are still several parameters such as the presence of vegetation above the soil, the swell, and shrink of soil that need to be investigated more in depth and the most important factor is the hydraulic properties of soil and its variability across landscape. This variable is more complicated, and there are many steps to reach an approximately correct value. In the current study, existing values were applied from previous work done around the same study site, while some other values were obtained from field and laboratory experiments. However, it is known that soils in agro-ecosystems tend to experience large changes in some of these properties and, therefore, are difficult to capture. This factor needs to be taken into account when implementing this model over unnatural setting systems.
Finally, based on current numerical model outputs and field experimental observations allowed identification new challenges in northern agro-ecosystems. Improved representation of soil dynamics is necessary to improve fidelity in the simulations, and also there is a need to establish better strategies to compare single-point numerical modeling to scale-dependent micrometeorological observations. In addition, a large deviation in simulated soil profiles and heat exchanges reveals the highly heterogeneous nature of an aerodynamically simple terrain considered in terms of atmospheric observations.