Modeling Soil Water Dynamics and Pasture Growth in the Montado Ecosystem Using MOHID Land

: The southern Iberian Peninsula is characterized by evergreen oak woodlands (locally known as montado ), which constitute an important savanna-type agro-silvo-pastoral ecosystem. This ecosystem is facing a progressive decline for several reasons, with the foremost being overgrazing. Better management tools are necessary to accurately quantify the systems’ carrying capacity and the sustainable stocking rates that prevent land degradation. The purpose of this study was to determine whether the MOHID-Land model could adequately simulate soil water dynamics and pasture growth in the montado ecosystem. The study area was located in the Alentejo region of southern Portugal. The model successfully simulated soil water contents and aboveground biomass during the 2010–2011 and 2011–2012 growing seasons, producing acceptable errors of the estimates (0.015 ≤ RMSE ≤ 0.026 cm 3 cm − 3 ; 279 ≤ RMSE ≤ 1286.5 kg ha − 1 ), and relatively high modeling efﬁciencies (0.481 ≤ EF ≤ 0.882). The model was further used to simulate the same variables for a longer period (1979/2009 seasons), to account for the effect of climate variability on model estimates. Water balance and dry biomass estimates were found to be signiﬁcantly different between rainfed and irrigated pastures, as well as between the ten driest and ten wettest seasons, with the model responding well to climate variability. The results showed the potential of using the MOHID-Land model for improving pasture management in the montado ecosystem.


Introduction
The southern region of the Iberian Peninsula is characterized by a savanna-type agro-silvo-pastoral ecosystem, known as montado in Portugal (hereafter adopted) and dehesa in Spain [1,2]. The montado system consists of an open formation of cork (Quercus suber L.) and holm oak (Quercus ilex rotundifolia L.) trees, presenting high levels of spatial variability in terms of density, combined with fallow land or pastures, which can be natural, improved, or cultivated [3,4]. The structural diversity of the system combines with the good overall habitat connectivity over its area of distribution, allowing animals to move around their territories, find mates, hunt, forage, and reproduce, resulting in high levels of biodiversity. As a result, the montado is considered as a High Natural Value system [2], providing various ecosystem services (ES) that are perceived by farmers, stakeholders, and society in general as being important for human welfare [5][6][7].
In the last decades, the implementation of the European Common Agricultural Policy has led to an intensification of grazing, which has caused additional stress to the system, promoting the spread of tree diseases and preventing the natural regeneration of the ecosystem [3,4,8]. The most significant transformations have been (i) the increase of the stocking (or livestock) pressure (i.e., the replacement of sheep by cattle, the replacement of light indigenous breeds of cattle by heavier breeds, and the increase of livestock units per area); and (ii) the use of heavy machinery for shrub control [7]. These changes have exhausted the natural pastures, decreased tree regeneration, and degraded water and soil quality [3,7]. Therefore, there is the need for the development of policy instruments that consider the specificities of the montado silvopastoral system; otherwise, its future will be severely threatened in the short term [5,7,9].
One major challenge for improving system management is understanding soil water dynamics in the pasture component of the montado ecosystem, allowing for its improved management regarding grazing, and hence its capacity to withstand increased stocking pressure. Mediterranean grassland species have evolved adaptive strategies to avoid or endure severe summer droughts, adjusting their lifecycle to seasonal water availability [10]. The emergence of plants depends on the timing of the first autumn precipitations, and their active period is in winter and early spring, with senescence in May [11][12][13][14]. In montado systems, the understory plays an important role in the functioning of the ecosystem [12], having a variable productivity that responds to the temperature and precipitation (P) regimes. The relationship between soil water and biomass also influences whether these ecosystems act as CO 2 source or sink, with Jongen et al. [12] suggesting that in dry years grasslands may act as a source of CO 2 , while in more humid years they may act as a sink. With climate change scenarios suggesting a decrease in annual precipitation in mid-latitude temperate regions, accompanied by increasing temporal variability of precipitation and drought risks [14][15][16], it is thus essential to develop tools for accurately quantifying the soil water balance and net primary production, and improving system management. The additional variability in precipitation also introduces the possibility of an additional management alternative: the use of of irrigated pastures. These both allow a smoothing of the effects of precipitation variability in rainfed pastures and provide a buffer of forage production, to release grazing pressure from rainfed pastures when that production is inadequately high.
Today, several process-based models, such as the BIOME-BGC [17], STICS [18], SPACSYS [19], EPIC [20], ARMOSA [21], and PaSim [22] models, are available for simulating soil water dynamics and vegetation growth in grasslands. Most of these models provide an integrated perspective by taking into account the various vegetation, soil, and weather processes regulating energy and matter exchanges in agroecosystems. Some of these processes though, as in the cases of the EPIC [20] and PaSim [22] models, are described in a simplistic way. For example, soil water dynamics is usually defined based on the notion of field capacity and the wilting point, with soil water storage capacity being computed from atmospheric demands, while other outputs, such as capillary rise and percolation, are computed empirically. Other process-based models, such as SPACSYS [19], and ARMOSA [21], make use of the Richards equation for computing soil water dynamics, thereby providing a more physically-based approach for predicting soil water contents and fluxes. Another example is the MOHID-Land model [23], which is a distributed model based on primitive equations for surface runoff and flow in both the vadose zone and the aquifer. Flow in porous media is based on the Richards equation, while infiltration can be simulated using Darcy's equation and surface pressure, due to the surface runoff water column or by using empirical algorithms (e.g., Green and Ampt). The MOHID-Land model further includes a modified version of the EPIC model [20] for simulating crop development and biomass growth, which can potentially be used for estimating the stocking rates in the montado ecosystem and prevent soil degradation.
Hence, the objectives of this study were (i) to calibrate/validate the process-based MOHID-Land model [23] for estimating soil water dynamics and dry biomass growth in a rainfed and irrigated pasture during the 2010-2011 and 2011-2012 seasons; (ii) to estimate pasture irrigation (I) needs, using the calibrated model for the 1979-2009 seasons (30 years); and (iii) to compare the soil water balance and dry biomass estimates in rainfed and irrigated pastures during both wet and dry seasons, providing an insight on the influence of climate variability on model performance. The hypotheses addressed in this study were (i) that the MOHID-Land model could accurately estimate the soil water balance and aboveground biomass growth in pastures under different management regimes; (ii) that the model could be used for irrigation scheduling; and (iii) that climate variability had a significant effect on model estimates. The results of this study will help to support the decision-making process regarding ecosystem management, by better accounting for the stocking rates and reducing the adverse impact of grazing in the montado ecosystem.

Field Site Description and Data
The field data was collected at Herdade da Machoqueira do Grou, located near Coruche, southern Portugal (39 • 08 16" N, 08 • 20 03" W). The climate in the region is dry and sub-humid, with mild winters and hot, dry summers. The mean annual temperature is 15.9 • C, with the mean daily temperatures at the coolest (December) and warmest (August) months being 10.4 and 23.8 • C, respectively. The mean annual precipitation is 680 ± 210 mm, with 87% falling between October and May (over the 1955-2007 period). The weather variables (rainfall (mm), average temperature ( • C), wind speed (m s −1 ), relative humidity (%), and solar radiation (W m −2 )) were provided for the study area by the MM5 mesoscale model (http://meteo.tecnico.ulisboa.pt), forced by the initial conditions from the NCEP (National Centers for Environmental Predition) Climate Forecast System Reanalysis, at a spatial resolution of 9 km.
The MOHID-Land model was calibrated/validated using data from Jongen et al. [13,14]. These authors carried out a rainfall manipulation experiment to study the resilience of montado understory to precipitation variability between October 2010 and June 2012. Eight rainfall manipulation shelters were constructed, each one covering an area of 6 m × 5 m (30 m 2 ), within a fenced area of about 3500 m 2 . The shelters were covered with a clear, 0.2 mm, UV-transparent polyethylene greenhouse film. More details related to the shelter design can be found in [13].
During the 2010-2011 and 2011-2012 growing seasons, the pasture within four manipulation shelters was irrigated every three weeks, receiving 40 mm per event, with this quantity being defined based on historical precipitation data . Hereafter, these plots will be referred to as "irrigated pasture". Water was applied to the experimental plots using an irrigation system, consisting of four 18-VAN rotary sprinklers with 90 • arc nozzles (Rainbird, Azusa, CA, USA), one in each corner of the experimental plots [14]. In addition, four non-sheltered "control" plots receiving natural precipitation were defined, referred to hereafater as rainfed pasture. All plots received an equal water input during the germination and seedling establishment (until the middle of November).
The soil volumetric water content was continuously measured (every 10 min) in the middle of each experimental plot with an EC-5 soil moisture sensor (Decagon Devices, Pullman, WA, USA). In the irrigated and rainfed plots, measurements were taken at 10 and 5 cm depth, respectively. Aboveground (ABG) dry biomass was measured four times per season. Two 30 cm × 30 cm quadrats were defined in each experimental plot. All plants inside were harvested, oven-dried at 60 • C for 72 h, and weighed. The maximum leaf area index (LAI max , m 2 m −2 ) was measured using a ceptometer (AccuPAR model LP-80, Decagon Devices, Pullman, USA), allowing an indirect determination of LAI by measuring the fraction of intercepted, photosynthetically-active radiation of the canopy. Maximum root depth was estimated at the end of the growing season (12 June). In each of the experimental plots, three soil cores of 8 cm diameter and a depth of up to 20 cm were taken. Roots were washed out, with subsequent analysis of root length for determination of specific root length (SRL) using WinRhizo software (Regents Instruments Inc., Québec, QC, Canada), then oven dried at 60 • C for 72 h, and weighed [14].

Soil Water Dynamics
MOHID-Land is an open-source, physically based, distributed model, which uses a finite-volume approach based on mass and momentum balance equations [23].
The variable-saturated one-dimensional water flow is described using the Richards equation: where θ is the volumetric soil water content (L 3 L −3 ), t is time (T), z is the vertical space coordinate (L), h is the soil water pressure head (L), K is the hydraulic conductivity (L T −1 ), and S is the sink term accounting for water uptake by plant roots (L 3 L −3 T −1 ). The unsaturated soil hydraulic properties are described using the van Genuchten-Mualem functional relationships [28]: where S e is the effective saturation (L 3 L −3 ), θ r and θ s denote the residual and saturated water contents (L 3 L −3 ), respectively, K s is the saturated hydraulic conductivity (L T −1 ), α (L −1 ) and η (-) are empirical shape parameters, m = 1 − 1/η, and is a pore connectivity/tortuosity parameter (-). The sink term (S) accounting for root water uptake in Equation (1) is computed using the macroscopic approach introduced by Feddes et al. [29]. In this approach, potential transpiration (T p ) is linearly distributed over the root zone, creating the function T p (z), and may be diminished by the presence of depth-varying root zone stressors, namely water stress [29][30][31]. For that, crop evapotranspiration rates (ET c , L T −1 ) are computed with the FAO Penman-Monteith method, using the single crop coefficient approach [32], and partitioned into potential soil evaporation (E p , L T −1 ) and potential crop transpiration (T p , L T −1 ) as a function of LAI [33]: where λ is the extinction coefficient of radiation attenuation within the canopy (-). Then, the actual transpiration rate (T a , L T −1 ) is obtained by integrating T p (z) over the root domain, and by limiting it using the piecewise linear model proposed by Feddes et al. [29]. This approach considers that the water uptake is at the potential rate when the pressure head is between h 2 and h 3 , drops off linearly when h > h 2 or h < h 3 , and becomes zero when h < h 4 or h > h 1 (subscripts 1 to 4 denote for different threshold pressure heads). E p is limited by a pressure head threshold value, in order to obtain the actual soil evaporation rate (E a , L T −1 ) [34].

Plant Growth
The MOHID-Land model includes a modified version of the EPIC model [20,35] for simulating crop growth. This model is based on the heat unit theory, which considers that all heat above the base temperature will accelerate crop growth and development: where PHU is the total heat units required for plant maturity ( • C), HU is the number of heat units accumulated on day i ( • C), i = 1 corresponds to the emergence date (-), mat is the number of days required for plant maturity (-), T av is the mean daily temperature ( • C), and T base is the minimum temperature for plant growth ( • C). Plant growth stages (initial, plant development, mid-season, and late season stages) are thus defined as a fraction of PHU. Crop growth is modeled by simulated light interception, the conversion of intercepted light into biomass, and LAI development [35,36]. The change in total biomass is calculated from the solar radiation intercepted by the crop leaf area, using Beer's law [37]: where ∆BM act,i and ∆BM i are the actual and potential increase in total plant biomass on day i (kg ha −1 ), RUE is the radiation-use efficiency of the plant ((kg ha −1 ) (MJ m −2 ) −1 ), PAR day,i is the daily incident of photosynthetically-active radiation (MJ m −2 ), λ is again the light extinction coefficient (-), and γ i is the daily plant growth factor (0-1), which accounts for water and temperature stresses. RUE is estimated using the approach proposed by Stockle et al. [38], while γ i is computed as Neitsch et al. [35]: where w strs,i and t strs,i are the water and temperature stresses for a given day i (-), respectively. Water stress (w strs,i ) is calculated as in Neitsch et al. [35]: where T a cum and T p cum are the cumulative T a and T p values on day i − 1 (note that, by construction, T a < T p , so water stress is bounded between 0 and 1). The effect of water stress on crop growth is thus only reflected on the following day. Temperature stress is computed as temperature divergence from the optimal (T opt , • C) using a sigmoidal function, influencing crop development similarly to water stress [35]: where T av is the mean air temperature for a given day ( • C), and T base is the plant's minimum temperature for growth ( • C), respectively. Leaf area index is computed as a function of heat units and plant stress [35]. During early stages (initial and plant development stages), LAI increment on a given day is a function of the fraction of the plant's maximum LAI (LAI max , m 2 m −2 ) that needs to be reached during those stages (fr LAImax,ini ), as well as crop stress. During the mid-season stage, LAI is assumed to be constant, while during the late-season stage LAI declines as a function of LAI max , PHU, and crop stress. Root depth is also computed as a function of heat units [35], while root biomass is assumed to decrease from 0.4 of the total biomass at emergence to 0.2 at maturity [39]. A more detailed description of MOHID-Land's crop module governing equations can be found in Ramos et al. [36].

Model Setup, Calibration, and Validation
In the MOHID-Land model, the simulation period covered the 2010-2011 and 2011-2012 growing seasons, being set from emergence to senescence each year. The soil profile was specified with 4 m depth and divided into three soil layers ( Table 1). The soil domain was represented using an Arakawa C grid type [40], defined by one vertical column discretized into 11 grid cells, 1 m wide, 1 m long, and with variable thickness (0.05 m on the top to 2.5 m at the bottom).
The upper boundary condition was determined by the actual evaporation and transpiration rates, as well as the irrigation and precipitation fluxes. ET c values were computed by multiplying the daily reference evapotranspiration (ET 0 ) values with crop coefficients (K c ) for the initial (0.3), mid-season (0.7), and late season (0.7) stages. These values were taken from Allen et al. [32], and correspond to standard K c values for pasture in the Mediterranean region. The K c value for the initial stage was then adjusted for the frequency of the rainfall events and average infiltration depths, while the K c values for mid-season and late season were adjusted for local climate conditionsm taking into consideration plant height, wind speed, and minimum relative humidity averages for the period under consideration [32]. Root water uptake reductions were computed by considering the following parameters: h 1 = −10, h 2 = −25, h 3 = −200 to −800, h 4 = −8000 cm [41].
Model calibration and validation were carried out during the 2010-2011 and 2011-2012 seasons, respectively, with procedures following Ramos et al. [36]. The soil hydraulic (Table 2) and the crop parameters (Table 3) were calibrated for both rainfed and irrigated plots (during the 2010-2011 growing season). A trial-and-error procedure was adopted. First, the parameters of the van Genuchten-Mualem equations (Equations (2) and (3)) were defined according to the average values proposed by Carsel and Parish [42] for each soil texture class. The soil hydraulic parameters θ s , α, η, and K s were then modified to reduce the deviation between the observed and the simulated soil water contents. The parameter θ r was not modified, as it did not significantly influence soil moisture simulations. The parameter was set to 0.5, as suggested by Mualem [43]. The crop parameters were calibrated by optimizing the default values of the EPIC model for pasture [35] until the deviation between the simulated and the observed ABG dry biomass was minimized. The values of the LAI max and maximum root depth were taken directly from Jongen et al. [13]. The calibrated model parameters were then validated during the 2011-2012 growing season and the statistical indicators (given below) computed. θ r , residual water content; θ s , saturated water content; α and η, empirical shape parameters; , pore connectivity/tortuosity parameter; K s , saturated hydraulic conductivity.
The goodness-of-fit indicators adopted for comparing MOHID-Land model simulations with the observed values of θ and ABG dry biomass were the mean error (ME), the root mean square error (RMSE), the normalized RMSE (NRMSE), and the model efficiency (EF). ME values close to zero indicate no bias. RMSE and NRMSE values close to zero indicate small estimation errors and good model predictions [44][45][46]. EF values close to one indicate that the residuals' variance is much smaller than the observed data variance, hence the model predictions are good. On the contrary, when EF is very close to 0 or negative, there is no gain in using the model [47]. Table 3. Optimized parameters of the crop growth model.

Crop Parameter Irrigated Plot Rainfed Plot
Optimal temperature for plant growth, T opt (

Simulation Scenarios
After model calibration/validation, the MOHID-Land model was used to compute the water balance and dry biomass yields in rainfed and irrigated pastures, hypothetically grown in the studied area during the 1979-2009 seasons (30 years). Climatic data was also provided by the same numerical mesoscale MM5 used earlier (http://meteo.tecnico.ulisboa.pt), and expresses the typical climatic variability found in the Mediterranean region (Figure 1), with annual precipitation amounting to between 68 and 805 mm.
In irrigated pastures, irrigation needs were further computed with a system-dependent boundary condition that automatically triggered irrigation when a certain threshold pressure head (h t ) was obtained in different grid cells of the root zone domain. Irrigation then ceased after a second target pressure head (h 0 ) was obtained in the same grid cells. This system-dependent boundary condition further included a series of constraints that prevented the application of meaningless irrigation amounts and countless irrigation events, namely a minimum irrigation pulse (I min ), a maximum irrigation pulse (I max ), and a minimum irrigation interval (I int ). The model was thus capable of automatizing irrigation whenever soil pressure heads dropped below h t in different cells of the root zone domain, supplying them sufficient water to reach h 0 in those same cells based on a pre-defined irrigation strategy [36]. Simulations for irrigation pastures were thus run with the following settings: h t = −800 cm, corresponding to h 3 in the Feddes et al. [29] model (i.e., no water stress was allowed); h 0 = −200 cm, corresponding to field capacity; I min = 5 mm; I max = 20 mm; and I int = 1 day. Results of the soil water balance components (E a , T p , and T a ), LAI, and dry biomass yields in rainfed and irrigated pastures were then compared for the 30-year period (1979-2019) and for the 10 driest and 10 wettest years ( Figure 1). Averages were compared using a paired, double-tailed Student t-test, in which a p value of <0.05 was considered significant. In the irrigated plots, θ showed a fast increase after rain or irrigation events, then decreased gradually due to redistribution, root water uptake, and soil evaporation. In the rainfed plots, soil water dynamics showed a similar behavior as in the irrigated plots. However, while θ was kept close to field capacity in the 2010-2011 season due to rainfall, in the 2011-2012 season θ was much lower, mainly due to the extended drought period that lasted from December 2011 to March 2012. Table 4  0.024 ≤ NRMSE ≤ 0.047), while the model efficiency was high (0.481 ≤ EF ≤ 0.863). According to Ramos et al. [48], Dabach et al. [49], and Sándor et al. [50], the deviation between measurements and model predictions may be attributed to several reasons, including model inputs, model structure, and field measurement errors, which can be considered also in this study. Nonetheless, the goodness-of-fit indicators were within the range of values reported in the literature for water content simulations using different process-based models [36,[51][52][53][54].

Pasture Growth
Pasture growth parameters were calibrated/validated for the irrigated and rainfed plots separately, based on the observed dynamics under the different regimes. However, only two parameters differed between the irrigation and rainfed regimes (Table 3): the fraction of PHU to reach the end of the canopy development stage (fr PHU,dev ), and the fraction of LAI max at the end of the canopy development stage (fr LAImax,dev ). These parameters influenced the time needed for pasture to reach the mid-season stage, with irrigated pasture reaching it faster due to higher water availability, and also presenting a higher LAI value at the beginning of that stage (i.e., higher transpiration rates).
Regarding the aboveground dry biomass, the MOHID-Land simulations were in agreement with the measured values ( Figure 3; Table 4). For the calibration period, the RMSE values were 1286.5 and 1125.5 kg ha −1 in the irrigated and rainfed plots, respectively. The NRMSE values were 0.210 and 0.372, and model efficiency were 0.584 and 0.869. For the validation period, RMSE, NRMSE, and EF values were within the same order of magnitude of those values. The general agreement between the simulated and measured data suggested that the calibrated parameters resulted in an adequate estimate of θ and aboveground dry biomass for both the calibration and validation seasons. Thus, it can be concluded that the model performance was generally good, and the model was properly calibrated for simulating soil water dynamics and crop growth in Herdade da Machoqueira do Grou. Table 5  . From that date on, plots were sheltered from rainfall, and an additional 370-454 mm were applied as irrigation [13,14] (Figure 2). Figure 4 shows the evolution of the T p , T a , and E a values in the irrigated and rainfed plots during the 2010-2011 and 2011-2012 seasons. T a varied between 143-166 mm in the irrigated plots, and between 56-65 mm in the rainfed plots (Table 5). These values corresponded to a water stress (given by the T a /T p ratio) between 0.86-0.95 in the irrigated plots, and between 0.60-0.94 in the rainfed plots.

Soil Water Balance
The largest water stress was obviously registered in the rainfed plot during the 2011-2012 season due to the extended drought period observed between December 2011 and March 2012. The E a values varied between 128-155 mm in the irrigated plots, and between 182-198 mm in the rainfed plots. E a values were thus notoriously higher than T a values in the rainfed plots. As discussed by Allen et al. [55], when revising the various methods for field ET estimation, combined approaches of accurate θ observations and water balance simulation modeling provide for appropriate accuracy in ET estimates. Nonetheless, those results are further consistent with Huxman et al. [56] and Kurc and Small [57], who referred to the fact that in semi-arid ecosystems E a may account for more than half of ET due to the near-surface source of E (the top 20 cm), and because small rainstorms only wet the top few cm of soil. Also, Graham et al. [58] reported pasture losing more water to soil evaporation when compared to ryegrass.  As the rainfall manipulation experiments were carried out between autumn and spring, most of the water applied was lost through deep percolation (226-362 mm), which was to be expected (Table 5), providing natural recharge to the existing aquifers.

Dry Biomass and Water Balance Estimates in Wet and Dry Seasons
Model estimates of aboveground dry biomass ranged from 1331 to 2681 kg ha −1 and from 2778 to 5092 kg ha −1 in rainfed and irrigated pasture, respectively, for the period 1979-2009 ( Figure 5). Mean differences between the two pasture regimes were found to be statistically different. Likewise, model estimates of the aboveground dry biomass for the 10 driest and the 10 wettest seasons were found to be statistically different. Mean values reached 2055 (1340-2469) and 4395 (3675-5092) kg ha −1 in the rainfed and irrigated regimes, respectively, during the 10 driest seasons, and 2132 (1476-2681) and 3802 (2778-4653) kg ha −1 in the same regimes during the 10 wettest seasons ( Figure 5). While mean values in the rainfed regime were relatively close, dry biomass estimates during the wet seasons were more regular, with the first (2067 kg ha −1 ) and third (2293 kg ha −1 ) quartiles being relatively close when compared with the dry seasons ( Figure 5). Note that simulations of dry biomass values started to differ between pasture regimes during spring, when water became a limiting factor (Figure 3). In these simulations, and based on the field data used for calibrating the model, irrigated pasture was modelled only until June/July, when plants reached their PHU. Therefore, pastures with longer life cycles-namely permanently irrigated pastures-will reach higher dry biomass values, since they are not limited by temperature and water stresses during spring and summer months. Irrigation was thus a determining factor for pasture's dry biomass increase. Irrigation needs varied between 20 and 360 mm, being higher during the dry seasons (220-360 mm) than during the wet seasons (20-200 mm). Irrigation needs were obviously dependent of rainfall (R 2 = 0.513), with the model trying to maintain soil water contents within the predefined threshold (h t ) and target (h 0 ) pressure heads. Figure 6 gives, as an example, the irrigation scheduling estimated by MOHID-Land for the 1980-1981 (P = 220 mm; I = 360 mm), 1995-1996 (P = 805 mm; I = 140 mm), and 2002-2003 (P = 592 mm; I = 20 mm) seasons. Results showed that rainfall amount and distribution had a notorious influence when computing the soil-water balance and irrigation scheduling with MOHID-Land. T p values were also found to be statistically different between rainfed and irrigated pastures. The difference was justified by the LAI values found in each pasture regime ( Figure 5), these also being statistically different. As explained earlier, LAI (Equations (4) and (5)) had direct influence on the partition of ET c values into E p and T p . Similarly, T a values ranged between 32-73 mm and 128-266 mm in rainfed and irrigated pastures (Figure 7), respectively, indicating a water stress that varied between 0.71-5.88% and 0.47-1.97% in the same plots. Water stress was generally higher during the 10 driest seasons, ranging from 1.42 to 5.88% in rainfed pasture and from 0.54-1.97% in irrigated pasture. During the 10 wettest seasons, these values tended to decrease (0.71-2.44% in rainfed pasture, and 0.55-0.99% in irrigated pasture). As in model calibration/validation, in rainfed pasture, E a values were always higher than T a values (83-256 mm). In irrigated pasture, E a values (133-293 mm) tended to be more similar to T a values (128-266 mm) (Figure 7), as plants' canopy considerably increased with the temperature increase during the final stages of pasture development (spring), covering the soil and decreasing the fraction of the soil surface exposed to radiation. During the wet seasons, E a values estimated for rainfed (118-156 mm) and irrigated (133-293 mm) pasture were more similar to one another than during the 10 driest seasons when the former decreased (83-185 mm) and the latter increased (172-252 mm).

Future Research Needs Direction
Although the MOHID-Land model results were relatively good, the modeling approach followed in this study was relatively simple, presenting still many limitations that need to be considered in future applications while addressing the complexity of the montado ecosystem.
Firstly, there is the need to simulate the carbon and nutrient fluxes, in order to better unveil uncertainties related to whether and when grasslands act as sinks or sources of ecosystem carbon, providing further insight into the effects of stocking pressures on soil quality. The current model formulation already includes a biochemical module to conduct those studies, and will be the focus of future simulations.
Secondly, there is the need to consider grazing removals during growing seasons, as the system should be seen as dynamic, and interactions with pasturing herds should always be taken into account in a proper management tool. Managing short-term variability in plant growth is considered to present a greater challenge in grazing systems than in cropping systems, which need to be tackled when developing a proper forecasting tool for accurately assessing daily/weekly pasture production [59]. However, for this study grazing was not considered, as the objective was to show the capability of the MOHID-Land model in simulating the main processes related to soil water dynamics and dry biomass growth in the montado ecosystem. Thirdly, the model currently considers grasslands as a single species, while in fact they should be treated as a mixture of plants with different strategic properties, as suggested by van der Molen et al. [60]. Thus, the model needs to be able to simulate different species within the same cell domain. This approach also needs to be adopted for modeling the rest of the montado ecosystem, namely the niche located directly beneath the tree canopy, and the inter-relationships between the understory and trees for soil water and nutrients [11]. Only by addressing this will it be possible to further understand the ecosystem decline observed in different regions of the southern Iberian Peninsula over the last few decades [6,7], namely the competition and synergies for water resources and nutrients.
Lastly, the MOHID-Land model is a process-based model that simulates variably saturated water flow using the Richards equation. This means that model simulations require reliable information on soil hydraulic properties for computing water flows, which are usually not easily available, as direct measurements are expensive, difficult, and time-consuming [61]. While there is the need to adopt indirect approaches, including pedotransfer functions [62], in order to provide the necessary inputs to process-based models, the approach followed in this study seems preferable to any of the existing simplifications, including those adopted by water-balance models. For example, Sándor et al. [50] compared nine simulation models for predicting θ in the topsoil and biomass production. These authors found that some of the main drivers and results of the grassland processes, particularly θ and yield, were not represented well by common grassland models, attributing deviations between observations and simulations to model formulations rather than structural uncertainties within the models.

Conclusions
The MOHID-Land model was able to successfully simulate soil water dynamics and pasture development in a plot located in southern Alentejo, Portugal. Model estimates of θ and aboveground dry biomass showed a general good agreement with field data collected during the 2010-2011 and 2011-2012 seasons. The RMSE, NRMSE, and model efficiency values were lower than 0.026 cm 3 cm −3 , lower than 0.047 cm 3 cm −3 , and higher than 0.481 cm 3 cm −3 for θ, respectively, and lower than 1286.5 kg ha −1 , lower than 0.450, and higher than 0.60 for the aboveground dry biomass. The MOHID-Land model was further shown to take into account climate variability when estimating the soil water balance and biomass growth in two different pasture regimes (rainfed and irrigated). Irrigation played a decisive role in the production of dry biomass, with irrigation amounts varying between 20 and 360 mm. As a result, actual transpiration values (128-266 mm) increased considerably when compared with the corresponding values estimated for rainfed pasture (32-73 mm), being closer to the actual evaporation values observed in both regimes (83-256 mm in rainfed pasture; 133-293 mm in irrigated pasture). However, the higher actual evaporation values show the inefficiency of the system, where most soil water is consumed through a non-beneficial use instead of being converted into biomass.
As precise quantification of water fluxes over the montado ecosystem is essential for accurately quantifying ecosystem carbon and assessing uncertainties related to the source or sink behavior of pastures, the MOHID-Land model can be considered as a valuable tool for farmers to take the stocking rates into account and reduce the adverse impact of grazing in the montado ecosystem.