Evapotranspiration Cycles in a High Latitude Agroecosystem: Potential Warming Role

As the acreages of agricultural lands increase, changes in surface energetics and evapotranspiration (ET) rates may arise consequently affecting regional climate regimes. The objective of this study was to evaluate summertime ET dynamics and surface energy processes in a subarctic agricultural farm in Interior Alaska. The study includes micrometeorological and hydrological data. Results covering the period from June to September 2012 and 2013 indicated consistent energy fractions: LE/R net (67%), G/R net (6%), H/R net (27%) where LE is latent heat flux, R net is the surface net radiation, G is ground heat flux and H is the sensible heat flux. Additionally actual surface evapotranspiration from potential evaporation was found to be in the range of 59 to 66%. After comparing these rates with those of most prominent high latitude ecosystems it is argued here that if agroecosystem in high latitudes become an emerging feature in the land-use, the regional surface energy balance will significantly shift in comparison to existing Arctic natural ecosystems.


Introduction
Recent warming in high-latitudes has significantly impacted Alaska's ecosystems [1][2]. These impacts have affected a broad spectrum of ecological, physical and societal systems of the Arctic [2][3][4][5][6][7][8][9]. In this context, agroecosystems and food related economic activities may be highly impacted by climate change over the next decades [10][11]. Therefore in order to meet future demands and conduct sustainable agriculture, considerable increase in food production must reduce the environmental impact [12].
Agroecosystems in Alaska currently represent only a small fraction of the entire landscape consisting mainly of boreal forests and tundra. Two current trends are supporting re-invigoration of Alaskan agriculture. First, desire for local food production and concerns about food security. Second, the combination of lengthening of the growing season, higher surface

Experimental Site
The field experiment was conducted at the Fairbanks Experiment Farm (FEF) on West Tanana Drive of the University of Alaska Fairbanks (UAF) Agricultural and Forestry Experiment Station (AFES), in Fairbanks, Alaska, USA (64°51 0 16.6@ N, 147°51 0 36.4 @ W, 150 m above sea level) (Fig 1). Experimental data were collected during the summer season from June to September of 2012 and 2013. The length of the growing season in the subarctic can be defined as the number of days between the last frost of spring and the first frost of fall. In this period of time the air temperature never drops below the freezing point [49]. Based on meteorological data covering the period 1906-2006, the length of the growing season in Interior of Alaska has increased over the last century by about 45% from 85 to 123 days [50]. The experimental site is characterized by an almost flat topography of the valley floor of the Chena and Tanana River basin. The area is sheltered on three sides from the northwest to the northeast by nearby hills rising to an elevation of 300-500 m, with another barrier about 250 km south the Alaska Range [51][52]. The site has three main vegetation types: woodland, grassland and crops, combined with bare land. The research area has a continental subarctic climate with long, cold winters and short, warm summers. Summer comprises the months of June, July, and August where air temperature average 15°C. On the other hand, winter months of November through March have average air temperatures of -18°C. Considering the central months of the two major seasons, the thirty-year (1981-2010) average air temperature in Fairbanks for July was 22°C with extreme temperatures rising above 32°C and decreasing down to -40°C in the month of December with continuous snow cover ground [53][54]. Precipitation is relatively low with the average annual accumulation for the period 1981-2010 about 263 mm, which mostly occurs in summer months of July and August [54]. In general, day length during the summer month rises up to 22 hours, which leads to a swing in temperature above 27°C for around 13 days. The range of frost free days, (i.e. air temperatures above 0°C) is approximately from 86 to 144 days with a median value of 115 days [54].
Approximately one hundred years ago the research site was part of the Alaskan boreal forest comprised mainly of Picea mariana specie and was cleared to be cultivated after 1906. In the past, large quantities of manure have been used as a supplement nutrient [27]. The land had long term tillage and crop residue management practices [55]. The site contains alluvial soil in the flood plains of the Tanana River. It is classified as a Tanana silt loam [53] with an approximate composition of 70% of silt, 8% of clay, and 22% of sand that has been conventionally farmed for about 80 years. It also contains a relatively high concentration of calcium carbonate and calcium sulfate at the soil surface [27]. A perched water table above the permafrost is about 8 m deep while the main water table is located about 20 m deep [27]. Crops were planted into the soil that had been summer fallowed the previous season. Some vegetable crops are usually grown with irrigation to improve and control crop growth allowing better use of the available plant nutrients [53]. As mentioned earlier this site was utilized by Braley [48] to estimate rates of ET from barley (Hordeum vulgare L) and rapeseed (Brassica rapa) fields during 1978 and 1979.
The research site for this study is considered a baseline for Interior Alaska agricultural research under the UAF AFES FEF and considered to be representative of floodplain Interior Alaska growing condition.

Instrumentation
Lysimeter setup. The experiment deployed several drainage lysimeters (S1 Fig) during the two growing seasons in 2012 and 2013 (June to September). Each lysimeter was built 62 cm long, 62 cm wide and 62 cm high. A lawn mix soil consisting mostly of sandy loam (66% sand), 29% silt, and 5% clay was added to each lysimeter up to 10 cm from the top of each lysimeter prior to the summer of 2012. The lysimeters were installed on a flat land area over a leveled horizontal plane. The bottom of each lysimeter had a 15 cm filter layer that consisted of stones, gravel and sand with a layer of geofabric above and beneath. The geofabric separates the soil from the filter layer. The filter layer kept the soil from spilling into the drainage and helped to drain water during heavy rains or irrigation events. A pipe collected the drainage from the lysimeter bottom.
In 2012, the experiment in lysimeter plots was conducted from 6 June to 16 September. A set of six drainage lysimeters were used to measure ET for lettuce (Lactuca sativa) after direct seeding on 8 June 2012. Sensors for soil volumetric moisture content (θ ly ) in the root zone were deployed at 15 and 30 cm depth in the lysimeters with a sample rate of 1 minute and record interval of 1 hour (S1 Table). Soil temperature at 15 and 30 cm depths was measured with sensors in the lysimeter (S1 Fig). The lysimeters were irrigated throughout the 2012 growing season. The observations were complemented by a multilevel 1, 2, 3 and 5 m meteorological observations, along with measurements of net radiation (R net ), turbulent velocities (u, v, w), and sonic temperature (T sonic ) operating from 16 July to 9 September 2012. Additionally, soil volumetric moisture content under barley, brome grass (Bromus inermis Leyss.), bare field (θ FEF ) and soil temperature (T soil ) at 15 cm depth were measured from 1 June to 30 September at the experiment site.
The intensive period of measurements during summer 2013 was from 14 June to 16 September, with two treatments in which ET was measured. The plot treatments (three replicates) were: (i) vegetated lysimeter and (ii) unvegetated lysimeter. A total of three soil moisture sensors were installed at 5, 10 and 20 cm depths in the vegetated lysimeters (θ ly ) and unvegetated lysimeters (θ unly ) treatments (S1 Fig). The soil temperature was measured at 5 and 10 cm depth in each lysimeter. Oak leaf lettuce (Lactuca sativa) at the five to six-leaf stage was transplanted into vegetated lysimeters on 1 June 2013. Irrigation was done throughout the 2013 growing season in both treatments by adding the same amount of water to each lysimeter. Irrigation amount ranged from 5.5 mm to 20 mm. Data from measurements in both lysimeter types were used three weeks after set up to allow the soil to settle. Observations during summer 2013 incorporated turbulent flux measurements derived from 3 m high sonic anemometer tower (operating from 7 July to 11 September 2013) and a large aperture scintillometer (LAS), which was operated from 7 July to 30 August. These observations were complemented by soil measurements (T soil ) in barley, brome grass and bare plots at 15 cm depth (operating from 1 June to 17 September) (S1 Table). All soil moisture and temperature profiles were recorded on dedicated data loggers (S2 Table).
Micrometeorological instrumentation. An eddy covariance (EC) instrument and meteorological station (Met station) were installed at the research site (Fig 1 and S1 Table). The EC instrument was placed on a tripod in the center of the farmland. The instrumentation consisted of a three-dimensional (3D) sonic anemometer (RMYoung 81000) mounted at a height of 3 m to measure the three turbulent components of the wind flow vector (u, w, v) with two temperature probes (T air ) mounted at 1 m and 3 m above the ground to determine air temperature (S2 Table). Data were collected at 20 Hz frequency and fluxes were calculated for a 30-min eddycovariance average period. With the aim to foster further studies a LAS was also installed on site (Fig 1). R net sensor was mounted at 3 m oriented to the south to avoid shade at all times. A barometer (P) was placed at the surface to determine the ambient air pressure. All data sensors were centralized in a single data logger (S2 Table). Additionally, two meteorological stations were mounted at 2 and 5 m above the ground surface to measure air temperature, relative humidity (RH), air pressure, wind speed (U), wind direction, and precipitation at 1-minute sampling rate. Data redundancy ensured a fairly continuous rate of data collection.
Pan Evaporation. A standard weather bureau Class A evaporation pan (PE) (122 cm diameter by 25 cm height), located 5 m away from the lysimeter plots, was used to measure manually (hook gage) and determine daily time series of potential evaporation (E P ). The water level in the pan was maintained within 7.5-12.5 cm of the lip. The evaporation pan is made of aluminum and rests on a wooden platform 12 cm above the ground over non-irrigated grass around the area. Daily Ep measurements were collected at 0800 AM AKST systematically every-day from 20 June to 5 September 2013 and corrected by wind observations atop the pan evaporation (S1 Table).

Surface Energy Balance
The surface energy balance is established based on Eq 1. As described previously the FEF sits on an almost flat surface terrain with no aerodynamic obstacles on the central section of the farm covering 1 km east-west and approximately 700 m north-south direction.
where R net is the surface net radiation flux (W m -2 ), G is conductive ground heat flux (W m -2 ), H is the sensible heat flux (W m -2 ), LE is the latent heat flux (W m -2 ) and Res is the residual closure component. In this case no storage term is considered since the vegetative canopy is very simple.
In Eq 1 the net radiation term R net , was measured directly. Ground heat flux (G) was calculated based on soil thermistors over different vegetation covers such as bromegrass field, barley field and fallow field. Soil temperature depths were 5 and 15 cm in summer 2012 experiment and 5, 15, 20 and 30 cm for the summer 2013. The conductive heat flux on the ground was calculated as follows: where T(z) is the soil temperature profile (°C) at specified depths z (cm) and k is the soil thermal conductivity. The k value in this study is treated as constant at 0.9 W (m°C) -1 [54]. The H component was measured based on meteorological data and compared to eddy covariance procedure. The latent heat term (LE) was estimated using ET method described in following subsections. Energy balance partitioning is used to determine the total available energy at the surface among the energy balance components by calculating the ratios LE/R net , H/R net and G/R net . These ratios indicated the relative magnitudes of LE, H and G in the surface energy balance. The ratio of H/LE flux is the Bowen ratio (β).

Energy balance closure
Based on independent measurements and determinations of R net , LE, H and G, the surface energy balance was established. Since Eq 1 combines radiative fluxes with turbulent fluxes averaged in space and time; still an energy closure was estimated characterizing the site in terms of the surface-atmosphere interactions. The closure fraction C F was therefore deduced based on Eq 3.

Estimation of evapotranspiration
Penman-Monteith. A number of approaches can be used to estimate ET based on energy balance measurements Penman-Monteith [56, PM hereafter], Priestley-Taylor [57] and Bowen ratio energy balance method [58][59]. Of these, the PM method is the more widely used in advanced ET models [60]. This method estimates ET on the basis of surface aerodynamic properties and physiological characteristics of vegetation. Variables used in the PM method are net radiation, soil heat flux, air temperature, relative humidity, wind speed, and environment-specific variables related to vegetation cover. The aerodynamic and physiological properties of the vegetation known as canopy resistance are the two important factors in the PM model. An approximation to actual evapotranspiration is the modified PM equation with the addition of the surface canopy resistance [61][62]. Several authors have shown the application of PM equation, including the canopy resistance, in a variety of environments. They have also tested this formulation including water stress, over several crops such as grass, lettuce, soybean, cattails, maize, tomato, wheat, and cotton [61][62][63][64][65]. Sensitivity of PM equation to different input data and parameters shows an effective dependence on the aerodynamic and canopy resistance introduced by considering the influence of vegetation type into the equation (see Eq 4) [66].
where ET is the latent heat flux of evapotranspiration (mm h -1 or mm day -1 ), λ is latent heat of vaporization (kJ kg -1 ), Δ is the slope of saturation vapor pressure versus temperature curve (kPa°C -1 ), R net is net radiation flux (W m -2 ), G is ground heat flux (W m -2 ), ρ a is the air density (kg m -3 ), C p is the air mass specific heat (kJ kg -1°C-1 ) at constant pressure, e s is the saturation vapor pressure at ambient air temperature (k Pa), e a is the actual vapor pressure of the air mass (k Pa), e se a is the vapor pressure deficit (VPD) (kPa), γ is the psychometric constant (kPa°C -1 ), r a is the aerodynamic resistance (s m -1 ), and r c is the bulk canopy resistance (s m -1 ). As indicated in Eq 4 the PM method requires information on net radiation, air temperature, air humidity, wind speed, and ground heat flux that can be obtained and deduced from meteorological and radiation observations. In this case, the vapor pressure deficit (VPD; e se a ) can be calculated as a function of measured air temperature and relative humidity using Eq 5 and Eq 6. where T is air temperature (°C) and RH is relative humidity in (%) The slope of the saturation vapor pressure (Δ) curve is also a function of temperature and can be calculated based on Eq 7. The psychometric constant (γ) in Eq (4) is a function of atmospheric pressure (which varies slightly over time and altitude) and is given by Eq 8: where C p is the specific heat at constant pressure, equal to 1.013 (MJ kg -1°C-1 ), λ is the latent heat of vaporization 2.45 (MJ kg -1 ), ε is the ratio of molecular weight of water vapor/dry air = 0.622, and P is the ambient air pressure (kPa). The aerodynamic resistance (r a ) under neutral conditions is calculated from Eq 9 following Allen et al. [67]: where r a is aerodynamic resistance (s m -1 ), z m (m) is height of the wind speed measurements, z h (m) is the height of temperature and humidity measurement, k is von Kármán constant (0.41), u z (m s -1 ) is wind speed measurement at z m , d (m) is zero plane displacement height of wind profile, z om (m) is roughness parameter for momentum, z oh (m) is roughness parameter for heat and water vapor. Reference values recommended in the literature are d = 2/3h c , where h c is crop height in meters; z om is 0.123h c , and z oh is 0.1 [68]. The stomatal resistance was obtained by measurements with a leaf porometer (S2 Table) in the field-scale to calculate canopy resistance. Following the procedure developed by Irmak et al. [65], we have randomly selected crops to determine the stomatal resistance in the field taking readings of leaves from the most representative vegetation patterns present in the farm (e.g., lettuce, barley, smooth bromegrass). Samples were taken during the midday interval from 1100 to 1400 AKST to determine a stable quantity representative of the central part of the day. Then, combining these data with the determination of the Leaf Area Index (LAI) measured by AccuPAR LP-80 ceptometer (S2 Table) the canopy resistance was determined. The variability range of canopy resistance was from 23 to 150 s m -1 with a median of 100 and standard deviation of 55 s m -1 . Thus, the median of 100 s m -1 was then applied into the PM equation to estimate the actual evapotranspiration.
Priestley-Taylor coefficient. Priestley and Taylor model (PT) [57] calculates potential ET based on the measurements of equilibrium evapotranspiration via an empirical coefficient α. This coefficient varies according to the surface and vegetation type. A constant value of 1.26 is generally used in landscapes where vegetation cover is almost complete and for saturated surface conditions [57,69]. The PT equation can be applied for unsaturated water surfaces provided α is adjusted to each condition [70]. The PT model has been shown to provide acceptable accuracy for predicting daily evaporation in Arctic ecosystems if the value of α is known [71]. The equation Eq 10 describes the PT approach. where α is an empirical coefficient relating actual evaporation to equilibrium evaporation; s is the slope of the saturation vapor pressure and air temperature curve (kPa°C -1 ); γ is the psychrometric constant (kPa°C -1 ); R net is net radiation (W m -2 ); and G is ground heat flux (W m -2 ). Mass balance approach. The ET estimated from lysimeters usually derives from applying the mass balance equation as a closed system as well as measurement of the soil water budget and some meteorological variables. The mass balance method is largely used in agriculture especially in crop productions that use irrigation input [67,[72][73][74]. The mass balance equation is indicated in Eq 11 [75][76].
where P is precipitation, I is irrigation, C r is capillary rise, ET is evapotranspiration that includes canopy interception or wet canopy evaporation and plant transpiration (i.e. dry canopy transpiration), D is drainage, R is runoff and, ΔS is the change in water storage (all terms expressed in mm) in both the unsaturated and saturated soil zones. In lysimeter systems, the runoff component R is not considered and the capillary rise C r is assumed to be negligible. The mass balance for the study can thus be expressed according to Eq 12: To calculate the lysimeter water storage the devices were divided into different layers for which measurements are available by depth, assuming that each soil moisture sensor was installed in a sampling depth layer within the lysimeter. According to Lewan and Jansson [77] on a similar setup, measurements at 5 cm depth were considered to represent 0-7.5 cm layer; 10 cm depth representing 7.5-15 cm layer, and 20 cm corresponding to 15 cm down to the bottom of the soil profile. Then, the value of the soil water storage was obtained per layer after integrating the volumetric soil water content in the specific depth. The total soil water storage was determined by the sum of the storage in each layer Eq 13 where S represents the storage (mm), θ the volumetric soil water content (m 3 m -3 ), and D the considered soil depth (mm). The change in soil moisture can be obtained by the change in soil moisture content over depth and time as indicated in Eq 14: where ΔS is the soil water storage, θ is the volumetric soil moisture content (m 3 m -3 ), t is the time, and z is depth (cm). Hence, the soil water storage variation profile was determined by the difference between the values of the soil moisture content obtained in the final and initial time of each considered period (daily or weekly), using Eq 15.
where ΔS is the soil water storage variation (mm), S f the final soil water storage (mm), and S i the initial soil water storage (mm). Two phases of crop developments (intermediate phase and maturity phase) were selected for the comparison between ET derived from mass balance and PM method. The period of 5 weeks after planting was identified as the intermediate phase, which started from 10 July to 23 July, 2013. The maturity phase was when the canopy is fully developed starting from 14-27 August 2013.

Meteorological and hydrological conditions
During the growing season, the average 30-min net radiation (R net ) in summer 2013 was slightly higher than 2012 (Table 1), with the daytime average of 156±122 Wm -2 in 2013 and 149±123 Wm -2 (Fig 2A) in 2012. Conductive ground heat fluxes (G) at the site were calculated according to Eq 2 and found to be mostly proportional to R net and following a diurnal cycle. On average, G in summer 2013 resulted to be similar to 2012 (Fig 2B). The mean air temperature was found to be on average 16.6°C in 2012 ranging from 0.2°C to 31.0°C and 18.2°C in 2013 with a variability range from -4.3°C and 34.9°C ( Table 2). The growing season 2013 included~5 days (44 hours) of negative air temperatures. However, the mean air temperatures in both years were higher than the 30-year average (Fig 3). The maximum air temperature reached to 31.0°C and 34.9°C during the month of June 2012 and 2013 respectively, while the normal (over 30 years) mean maximum was only 21.2°C (Table 2). This increase in maximum air temperature indicated a slightly warmer growing season in this high latitude agroecosystem. On the other hand, the minimum air temperature occurred in September with the lowest value of -4.3°C recorded in 2013.
The relative humidity (RH) of the experimental site averaged 69±19% for 2012 compared to 66±21% for 2013. High values of RH were correlated to lower air temperatures ( Fig 2D) as well as to increased precipitation events. Half-hourly mean vapor pressure deficit (VPD) varied during the growing season as depicted in Fig 2E. The mean midday VPD was 0.6±0.5 and 1.2 ±0.9 kPa with the maximum VPD of 3.4 and 4.3 kPa in June 2012 and 2013, respectively ( Table 1).
The precipitation field was found to be very variable and significantly different from long term averages. Collected values at the experimental site for both years (Fig 2F) resulted in much lower amounts than those from 30-year average 165.3 mm (Table 2). On the other hand, comparing side-by-side both summers it was found that August 2013 (56.2 mm) verified larger amounts than August 2012 (30.5 mm) while the normal monthly average precipitation is (47.7 mm). Similarly, the driest period in the past 30 years was verified to be June of every year (mean precipitation of 34.8 mm); however, June 2013, showed in the study area, a precipitation of 4.4 mm that was below the 30-year average. In contrast, the amount of precipitation in June 2012 (53 mm) was higher than the normal average of June. In addition, precipitation decreased about 6% during the summer season of 2012 when compared to the long-term 30 years mean  ET-Cycles in Northern Agroecosystems for this area [16]. A decreased rate of more than 28% was found in the summer 2013 resulting in abnormally dry conditions. In terms of T soil , there has not been significant differences in average T soil measured at 15 cm depth at the experiment site when comparing the growing seasons of 2012 and 2013 (Fig 4; Table 1). The soil volumetric water content (θ ly ) in the lysimeters in the layer 0-15 cm (2012) and in the 0-20 cm depth (2013) varied greatly over the growing season. The variability of θ ly depended on irrigation practice (i.e., irrigation quota and timing) and precipitation events. An average θ ly was 0.3851±0.0174 (6 June to 27 July 2012) and 0.3685±0.0245 m 3 m -3 (1 June to 6 September 2013) while an average θ FEF in the farm field with no irrigation was 0.1700±0.0197 m 3 m -3 average accounting for farm diversity of land surface type (e.g., crop land, grass land and bare land) ( Table 1).
The frequency distribution of surface wind direction and wind speed is shown during the period of study for the 2012 (

Surface Energy Balance
Energy closure at half-hour time-scale. The energy balance closure is defined as the ratio between the resulting turbulent fluxes manifested at the surface and the total energy available [80][81][82]. In this study, the energy balance closure was evaluated for the entire dataset 1,540 thirty-minute daytime intervals. There is a strong linear relationship between the sum of the 30 min average latent heat (LE) and sensible heat (H) plotted against the available energy (R net -G) for the summer 2013 growing season (Fig 6). A slope of 0.95 and an intercept of 10 W m -2 was obtained. These values indicate that on average the turbulent heat fluxes are slightly underestimated (by~5%) neglecting the storage term in the energy balance equation due to the short canopy across the farm landscape (i.e., less than 0.50 m on average; see Fig 1). Similar results were found by Li et al. [83] over maize farmland in Northwestern region of China, while Parent and Actil [84] obtained 0.79 in a farmlands at Saint-Ubalde, South-Eastern Canada. Moreover in grasslands, energy balance closure values generally ranged from 0.70 to 0.80 [80], 0.70 in prairie [85], 0.74 in olive orchard field [86], 0.77 in switchgrass field [87], and 0.85 in an alfalfa field [88]. However in terrestrial ecosystems, particularly including forest [82,[89][90][91][92][93] the energy balance closure was found to range from 0.50 to 0.96 due to the complexities of the canopy architecture.
Therefore, it is concluded that the surface and atmospheric flow conditions obtained in this study (Fig 1) established the energy balance closure that is highly reliable and useful for examining energy partitioning among all energy fluxes.
Energy balance and energy partitioning.  reached the peak value of 296 W m -2 around midday, basically following the time-variation of R net . Then, LE rapidly decreased to zero at 2100 AKST when transition in the atmospheric surface layer started. On the same day, H flux slowly rose from 0 to a value of 180 W m -2 , when R net peaked at 385 W m -2 , then H declined steadily to zero at 1900 AKST indicating the changes in stability conditions in the atmospheric surface layer. The midday fraction of available energy (R net -G) into H was about 37%. On the same day, the G flux was the smallest compared with the rest of fluxes and became positive at 1000 AKST. The peak magnitude of G was 38 W m -2 and occurred at roughly 1600 AKST. Later on G dropped below zero about 2200 AKST~two hours after R net turned negative. The calculated Bowen ratio (β) around middday was 0.6.
The monthly midday averages of the energy balance components (R net , LE, H, G) for two years were calculated, and these results are illustrated in Fig 8. The data were acquired in the period June to September of each growing season. The summer mean R net flux input for both years peaked in June (reaching~244 W m -2 in 2012 and 264 W m -2 in 2013) and dropped-off gradually in August below 31-35% and September below 56-64% from the seasonal maximum (Fig 8). The amounts of R net for both years were slightly different, with~20 W m -2 in 2012 being lower than in 2013. When the midday means were considered over the period of the   On the other hand, it was observed over the two years that the major portion of surface energy balance was attrributed to LE in the period June to September behaving similarly to that in tundra and wetlands in Arctic and Subarctic sites (Table 4). Nevertheless, LE/R net value is comparable to the lower range variability of that obtained from maize and soybean farmland, Nebraska, USA (0.6 to 0.9) [94] and are close to the value of 0.68 obtained in commercial farms near Flora city, Florida USA [95]. In addition, the value of H/R net ranged from 0.28 to 0.37 (mean 0.33) and 0.20 to 0.30 (mean 0.25) of R net in 2012 and 2013, respectively. In comparison to other farmlands the ratio of H/R net in this study was found slightly higher than in soybean and maize (0.2 to -0.2) according to experiments carried out in Nebraska [94] considering only periods of fully-developed canopies.

ET-Cycles in Northern Agroecosystems
At all instances the β in the experiment site during the summer 2012 and 2013 was found to be systematically less than unity and ranging from 0.3 to 0.9 with an average of 0.36 in 2012  First column represents the ecosystem types, second column is the location of measuring site, third to fifth columns are energy partitioning values for LE /R net , H /R net , G /R net (derived from daily midday flux averages), sixth column is VPD (kPa) for each ecosystem type, seventh column is the Bowen ratio (β), eigth column is the measuring method used for energy budget components measured, and nineth column is the reference for data. EC = Eddy covariance method. BREB = Bowen ratio-energy balance method.  (Table 3). On the other hand, G/R net ranged from 0.05 to 0.11 in 2012 (mean 0.06) and 0.01 to 0.08 (mean 0.07) in 2013 of R net . Large values of G are typically found in subarctic landscapes [95] with the exception of boreal forests [96]. Nevertheless, the ratios of H/ R net and G/R net were in the interval of observed values from wetlands and shrub tundra (H/R net near 28%, G/R net ranges from 6-12%) from the Western and Central Canadian Subarctic [97].
In terms of evaluating the general trends on ET associated with changes in vegetation, soil moisture, and meteorological parameters [69] the Priestley-Taylor coefficient (α) was calculated. High values of α are associated with a high-energy partition of LE/R net , while low value represents the opposite. Stewart and Rouse [98] found that the theoretical value of α = 1.26 is generally applied to saturated surfaces in high latitude. However, in the present study α = 0.91 was the average from two years while α ranging from 0.40 to 1.22 is reported for Arctic and boreal ecosystems [44]. The average value of 0.91 in this study is consistent with values reported by Eaton et al. [97] in upland tundra. This variable range of α depends on the specific ecosystems under consideration for example, Liljedahl et al. [99] reported mean midday α = 1.08 (offshore) and 0.95 (onshore) in Arctic coastal wetland.

Evapotranspiration from water balance equation
ET by mass balance from irrigated lysimeters versus ET by energy balance. ET based on energy balance was obtained using Penman-Monteith (ET PM ) approach. ET PM was computed based on available meteorological variables collected at the site (Fig 2), Eq 4 according to Monteith [100]. The cumulative ET from mass balance obtained by Eq 12 was applied to the irrigated vegetated (ET VL ) and unvegetated (ET UVL ) lysimeters. The measurements of precipitation (P), irrigation (I), drainage (D), and change in storage (ΔS) allowed estimating ET.
Comparing ET rates among different treatments on lysimeters, the results showed ET VL having from 5 to 25% larger cumulative ET compared to ET UVL . When considering ET PM as reference of a larger evaporative area, the ratios ET VL / ET PM and ET UVL / ET PM verified a lower fraction than 1 during the first week. While, for the rest of the experiment, it resulted in a ratio larger than or fairly close to 1 (Table 5). Similar results were obtained by Braley [48] based upon their study within irrigated lysimeter and non-irrigated lysimeter in 1979. On the other hand, the ratio of ET VL / ET PM was found to be slightly higher than ET UVL / ET PM (Table 5). Nevertheless, on average ET mass balance was mostly higher than ET energy balance due to additional water input from irrigation. The average ratios of ET VL / ET PM and ET UVL / ET PM were found to be 1.12 and 0.97 , respectively. However, the ratio of ET from mass balance to the measurement of pan evaporation averaged 0.59 and 0.66 for ET UVL and ET VL, respectively. First column is time period covered by measurement (i.e., 10-16 July, 2013), second column is a weekly accumulated ET in the vegetated lysimeter (ET VL ), third column is a weekly accumulated ET in unvegetated lysimeter (ET UVL ), fourth column is a weekly accumulated ET by energy balance (ET PM ) derived using the Penman Monteith equation, fifth column is the ratio of ET VL / ET PM , and sixth column is the ratio of ET UVL / ET PM . doi:10.1371/journal.pone.0137209.t005

ET-Cycles in Northern Agroecosystems
Additionally, the ET estimate from water mass balance approach provided higher rates than from energy balance approach, and this difference was accentuated as the vegetation fully developed (Table 5). However, ET from energy balance method can be used as a reference ET for an agroecosystem especially in the sparse vegetation landscape.
In order to compare and benchmark the hydrological rates in agricultural lands, Table 6 shows the annual and summer hydrological balance characteristics among various ecosystems in Arctic and Subarctic regions. Based on, the total precipitation of 65 mm and irrigation 41.2 mm during this study period, the ET VL was almost 97% while the ET UVL was approximately 88% of precipitation and irrigation. In contrast, lower percentages indicated in Imnavait Creek Basin in North Slope of Alaska reported that 50% of precipitation went through the ET process and only 36% was found in the Upper Kuparuk Alaskan watershed [101] In addition, 76% of precipitation was found to be evaporated from the permafrost in the boreal forest at Caribou-Poker Creek Watershed in Interior Alaska [102]. Nevertheless, other studies have also shown a lower ratio of precipitation being evaporated through ET process when compared to this present study [103][104][105][106][107][108][109] (Table 6).
Penman-Monteith Evapotranspiration (ETPM) and Pan Evaporation (EP). Potential ET was measured with a Class A evaporation pan (E P ). Daily pan evaporation determinations (E P ) were manually made at 0800 AKST and no later than 0815 AKST every day. The E P fraction is defined as the potential evaporation rate for a given location. E P in this study ranged from 0 to more than 8.57 mm per day under clear skies conditions with daily average of 3.44 ±2.15 mm per day. Because manual E P measurements were made at one given time everyday the temporal series of ET PM was then compiled for a similar time interval for daily estimates comparison. There were about 89 measured values of E P available and only 69 values were used for comparison with ET PM because of sensor malfunctioning (S1 Table). On the other hand, during the study period, relatively high rates of ET PM were recorded in the month of July, while a declining trend was shown in September (Fig 10). Daily values of ET PM ranged from less than 1 mm to more than 4 mm and the daily average was 2.27±1.40 mm per day. This average value is slightly higher than values obtained by Braley [48] for the same site. The results showed that overall daily values of E P exceeded ET PM . Regression analysis was used to relate ET PM to E P and a correlation (R 2 = 0.69) was found, while a poor correlation (R 2 = 0.38) between those values was documented in other environments Florida, USA [95].

Discussion
This study investigated the energy and water mass balance during 2012 and 2013 growing season at the UAF AFES FEF representative of high latitude agricultural lands in Interior Alaska. The summertime climatic characteristics at the site during both years were examined. An increase in the mean air temperature of +2.2°C and +3.7°C was observed in 2012 and 2013, respectively, when compared to the 30-year average of mean air temperature. Nevertheless, the mean air temperature regime during 2012 resulted within the normal range of variability for 30-year climatological data ( Table 2). It is worth mentioning, that the mean values during 2013 verified a temperature excursion larger than one standard deviation when compared with the 30-year climatological mean. This warmer mean air temperature observed in summer 2013 is consistent with recent results indicating an increase in air temperatures of 1.4°C for Interior Alaska during the last 100-year record [16]. On the other hand, summer precipitation for 2012 remained approximately within the range of the 30-year average (165.35 mm); while, precipitation for 2013 was~26 mm below the normal average. After statistical examination of time series of R net values it was concluded that no major differences between years were found ( Table 1 and Fig 8).
In terms of turbulent flux regimes, the sensible heat verified no major variations with an average of~80 W m -2 . However latent heat fluxes increased by~22 W m -2 during summer 2013. This slight positive trend could only be explained by the change in pre-season conditions that made the sub surface to be wetter through an extended snowmelt period than in the other year (2012). [79][80].
The value of energy balance closure C F found in the field experiment reached levels of~0.95 (Fig 6). This value is considered to be representative of energy closure in agricultural fields because often the topography characterizing such systems is close to ideal conditions (i.e., flat terrain covered by short grass). Similarly, we found this value to be in good agreement to closure levels 0.70 to 0.99 observed at Fluxnet sites including several agricultural lands [91]. Still, a small residual term was found~3%. This term is generally attributed to systematic methodological errors, systematic instrument bias, neglected energy sinks, and unrepresentativeness of the G term [81,91,93,[109][110][111][112][113]. After a careful revision of all terms intervening in the energy balance the residual term, can only arise from surface patches containing different crops (e.g., bare land, grassland and trees). Therefore, an evaluation of G was conducted over the mentioned surface patches and it was observed to vary 3-6% (Table 1). This attribution is in agreement with other reports [93,114] in which G was found to dominate the relative uncertainty on the surface energy balance closure reaching up to 20% in agriculture sites [115][116].
The energy partitioning of R net into H, LE and G is strongly influenced by changes in surface conditions such as dynamics of vegetation growth, changes in soil moisture and surface temperature affected by precipitation [117][118][119]. In particular the relationship between LE and soil moisture is complex, variable in space and time and verifies nonlinear relationships with the energy balance terms. For example, LE/R net was found practically similar in dry and wet soils constrained in the case by VPD<0.30 kPa in Arctic coastal wetland [99].
In the present study, the energy ratio of LE/R net was found to be systematically larger than the ratio for H/R net and G/R net consistently also with β < 1 (Table 3). Likewise, several ecosystems in the Arctic and Subarctic have larger LE /R net than H/R net [44,45,97,[120][121][122][123][124] (Table 4). Alternatively, we have found that in comparison to some Arctic ecosystems [44,97,99,[125][126] this ratio LE /R net is largest amongst the other fractions. However, we have to point out that this ratio is still on the lower range interval when compared to mid-latitude agricultural fields [127].
The monthly trends of energy fractions accounting for their seasonal evolution were observed to maximize around the middle of the summer to then trend negatively to the end of the season (Fig 9). This behavior is verified in the case of G/R net and LE/R net . However, the energy ratio H/R net is observed to fluctuate at the end of the season in 2012. These variations in H/R net are consistent to changes in the thermodynamics of the air mass as indicated in Fig 4. Similarly the energy fraction associated to LE /R net verifies a positive trend during the decaying phase of the season demonstrating an increasing response to monthly precipitation during August 2012 (Table 2).
With the aim to identify agricultural land energy fractions in the framework of natural ecosystems in high latitudes, Table 4, reports a comprehensive comparison among these systems across the panArctic. For Arctic and subarctic wetlands, LE/R net was reported to be larger than 0.57 according to the studies of Moore et al. [128], Harazono et al. [119], and Rouse [121]. In contrast, a lower value of LE/R net in Arctic coastal wetlands was observed. In this case, a different environmental forcing due to the presence of onshore wind constantly offsets the energy partitioning [99]. While on the other hand, the reported values in literature of H/R net and G/ R net are similar to the ones in the present study.
Furthermore the energy partitioning in tundra ecosystem verifies in comparison mostly a lower LE/R net ranging from 0.36 to 0.67 [44,91,120,122,129]. On the other hand, H/R net in the present study compares well with the lower range reported from the mentioned studies in the range 0.26 to 0.40. Finally the energy fractions obtained in this study correspond well with the results obtained in upland tundra ecosystems reported by Eugster et al. [44] in which LE was the dominant component of surface energy balance. Conversely, energy balance studies in Alaskan coniferous boreal forest (i.e, composed mainly of white and black spruce trees) have found that H dominated the energy balance [96,122,130) with the exception of the study of Nakai et al. [123] in the Poker Flat Research Range which indicated LE slightly dominant on a sparser canopy over discontinuous permafrost.
In terms of evaluating the ET rates, this study produced two different approaches based on mass balance (i.e., lysimeters based) and energy balance (i.e., micrometeorological based). These two approaches have definitively different spatial and temporal scales in terms of their environmental interactions [131] and therefore their ET rates were slightly different owing to the vegetation development and the spatial scale representation. In order to evaluate the potential for environmental interaction of agricultural land the lysimeter experiment was conducted based on irrigation practices over two treatments: vegetated and unvegetated. Overall this study found that ET VL was higher than ET UVL while ET UVL was similar across the season to ET PM (Table 5). However, it is important to note that if ET PM is taken as the reference, the ratios ET VL /ET PM and ET UVL /ET PM verified a lower fraction in the first week due to the development phase of the vegetation. On the other hand, the fraction ET UVL /ET VL represented the percentage of ET due to vegetation growth and interception ranging from 75 to 94% (Table 5).
To give a prospective impact of agricultural lands in the framework of high latitude environments, Table 6 brings similar data records around the pan-Arctic from ten sites which are compared to the findings of this study. Establishing the ratio ET/P allows the evaluation of the percentage of precipitation input that is sent back to the atmosphere through ET. Based on the synthesis of Arctic basins hydrology study by Kane and Yang the ET/P ratio is well-correlated to latitude in Arctic natural ecosystems and basically accounts for 36-75% of the mass balance (Table 6) [101][102][103][104][105][106][107][108][109]. Whereas, in this study this fraction ranged much higher from 87% to 97%; only comparable to the ones obtained by Thorne and Hawkins [106] calculating 80% return (Table 6). These differences in ET and ET/P ratios are due to availability of energy for fluxes at the surface [132][133], precipitation distribution and rate as well as topography [134], forest canopy interception capacity associated to tree species and leaf area index [135]. It is therefore concluded here that the fraction of ET returned to the atmosphere in agricultural lands represent a much larger fraction of what has been reported for boreal forest at approximately the same latitude [96] and also larger than the fractions obtained in Arctic tundra.

Conclusions
We found that the ET cycles represent a large portion of surface energy balance partitioning accounting for approximately 67% of the net radiation. The ratio of ET obtained by water mass balance related to the measured potential ET ranged from 0.59 to 0.66 for evapotranspiration rates based on unvegetated and vegetated lysimeters, respectively. Additionally, ET was responsible for removing 97% and 88% of the moisture added to the vegetated and non-vegetated lysimeters, respectively. Northern latitudes are characterized by diverse ecosystems where wetlands and tundra dominate Arctic regions, and boreal forest with coniferous and deciduous trees dominate Subarctic regions. This work puts in perspective and compares the surface energy fraction on agricultural lands in the context of boreal forests, Arctic wetlands and tundra (Tables 4 and 6). The results indicate that the energy fluxing regime in terms of ET/R net of agroecosystems in the subarctic exhibits similar characteristics to tundra in the Arctic; contrasting therefore with subarctic boreal forest. Moreover, differential fluxes may manifest between agricultural and boreal forest over short spatial scales forcing small-scale circulations creating an additional imbalance term in the energy budget. Therefore, this study indicates that the presence and further development of agroecosystems in northern latitudes may lead to significant changes of ET cycle during the growing season in comparison with natural existing ecosystems.
Consequently, replacing native ecosystems to promote agricultural development and economic activities may result in significant changes in the land-use and therefore in surface energy regimes and balance. Moreover, these changes can collectivelly upscale to shift seasonal magnitude (season length and timing) and temporal partitioning of regional fluxes introducing a positive feedbak to climate.
Finally, on the basis of a changing climate scenario manifested through increasing air temperatures, lengthening of growing season and changes in vegetation gradients in northern latitudes, expanding agricultural lands may lead to an increase of ET cycles (water vapor return to the atmosphere) that will propagate to larger atmospheric scales through nonlinear interactions characterizing the surface-atmosphere system.