A process-based numerical approach to estimate forest groundwater consumption in flatland petrocalcic soils

Flatlands host important economic activities as is the case of the Pampean region. It is one of the most important landscapes in Argentina, reaching an annual gross margin of about 146 U$D ha⁻¹ in 2001–2005 (Carreño et al. 2012). As in other regions of the world, natural grasslands are being replaced at increasing rates by forests, either by commercial plantations or tree invasion (Richardson 1998; Geary 2001; Herron et al. 2002; Jackson et al. 2002; Zalba & Villamil 2002; Paruelo et al. 2006; Jobbágy & Jackson 2007; Nosetto et al. 2008, 2012). This fact has gained importance in the last several years because of the large tree-products demand, the carbon sequestration market, and particularly by federal and provincial tax and financial incentives (Wright et al. 2000).

Changes in land use as transitions between grasslands and forests could alter the water budget and soluble salt fluxes (Jobbágy & Jackson 2004). Trees affect water fluxes and influence the direction and intensity of salt exchange between ecosystems and groundwater (Calder 1998; Jackson et al. 2001; Gyenge et al. 2002; Nosetto et al. 2008, 2012). In some conditions, it can also lead to an intense accumulation of salts in soils and aquifers (Jobbágy & Jackson 2003, 2004). At least three processes are suggested to be responsible for these chemical changes (Jobbágy & Jackson 2007): (1) the intrusion of deep salty water in the vadose zone due to tree consumption, (2) solute accumulation in the soil as they are excluded by root membranes, and (3) solute accumulation toward the soil surface as they are recycled by the ecosystem.

At a regional scale, a theoretical framework controlling the potential occurrence of soil salinization was proposed (e.g., annual climatic water balance, hydrological factors, etc.) (Nosetto et al. 2008). Thus, soil salinization is even more intense in regions with a precipitation deficit due to water requirements of tree plantations, and where groundwater use compensates this deficit (Nosetto et al. 2008). The occurrence of soil salinization by phreatophytes was observed in several regions around the world (e.g., Thorburn 1997; Heuperman 1999; Nosetto et al. 2007), even in the sub-humid Pampean region (e.g., Eucalyptus camaldulensis, Jobbágy & Jackson 2004, 2007; Nosetto et al. 2008). However, at a local scale, the potential occurrence of soil salinization is not straightforward. Even if the aforementioned conditions for soil salinization are reached, other factors as lithology (e.g., the presence of massive bedrock or cemented petrocalcic horizons) could inhibit soil salt accumulation (Jobbágy & Jackson 2003; Nosetto et al. 2008).

Water flow in flatland areas is dominantly vertical due to its negligible slopes that do not promote surface runoff (Varni & Usunoff 1999; Nosetto et al. 2007; Devito et al. 2012). Most of these systems around the world show a widespread presence of heterogeneities as petrocalcic horizons at different depths, thicknesses, and degrees of cementation (Duniway et al. 2010; Kuznetsova & Khokhlova 2015). Petrocalcic horizons could modify the soil water flow patterns and plant available water in afforested areas because they are traditionally considered as a barrier to root growth and water flow (e.g., Shreve & Mallery 1932; Pazos & Mestelan 2002; Dietrich et al. 2014).

Water table fluctuations are habitually used to quantify the forest groundwater consumption (e.g., Salama et al. 1994; Engel et al. 2005). This methodology is based on the direct use of diurnal fluctuations, or bore hydrograph separation approach overcomes the problem of extrapolating water use because changes in groundwater levels beneath a forest represent an integration of the water balance components in groundwater beneath the area. However, this approach is suitable to be applied in (1) those sites closer to the groundwater divide to minimize groundwater flow into the forest and (2) aquifers with low hydraulic conductivities (i.e., 0.01–0.1 m day⁻¹) to restrict the attenuation of groundwater fluctuations beyond the forest (Salama et al. 1994). However, the application of this methodology seems to be in appropriate in sites with (1) mean aquifer hydraulic conductivities higher than 5 m day⁻¹, (2) gradients <0.2% (e.g., see Varni & Usunoff 1999) or (3) low forest groundwater usages.

Process-based numerical modeling (i.e., water and vapor fluxes, and energy transport) in the soil–plant–atmosphere continuum appears to be a tool capable of capturing the coupled physical/biological processes controlling the flow processes in afforested flatland areas (e.g., Min3P model, Mayer et al. 2002; Bea et al. 2012). Thus, this work is devoted to present an alternative methodology based on process-based numerical modeling to quantify the forest groundwater consumption and plant available water in sites with the presence of highly cemented petrocalcic horizons. The proposed procedure was applied in monitored (i.e., water contents, temperature, sap flow, etc.) grassland/forest plots with these features located in the sub-humid Pampean region in Argentina. The obtained results could clearly be extrapolated to other sites with similar features. Since the salinization process is postulated as the main consequence of grassland afforestation, another objective is to assess the risk of occurrence in these environments.

The site of study is located inside the southern of Azul Creek watershed that lies within the Chaco-Pampean plain between 36°49′39″ and 37°21′6″ S, and between 60°9′10″ and 59°45′53″W (see Figure 1). This watershed shows average surface slopes of about 5%, 0.5–0.8%, and 0.2% in the upper, middle, and lower part of the hydrological basin, respectively (Varni & Usunoff 1999). Both grassland and afforested plots are placed in the upper part of the watershed.

The climate of this part of the Pampean region was classified as temperate, sub-humid/humid (Burgos & Vidal 1951). In contrast to other regions located at comparable latitudes in the Northern Hemisphere, it is characterized by an east-west moisture gradient and increasing continentality toward the northwest (Burgos 1968). Average temperatures in the northeast part of this region are 24 °C (January) and 10 °C (July) in summer and winter, respectively, whereas in the southwest zone, they are 20 °C and 7 °C, respectively (Prieto 1996). The sedimentary basin is filled with Late Pleistocene–Early Holocene loess and aeolian sands forming a continuous mantle across these sites (Zárate 2003). They present strongly cemented interbedded soil horizons mainly crust composed of carbonate minerals (0.7–1 m thick, see Figures 1(a) and 2; Zárate & Mehl 2010). They are covered by a mantle of very fine sand, silty loam and clay, as well as friable sand-clay loam sediments. They correspond to loessic deposits that constitute the parent material of a considerable part of the current soils; for this reason, they exhibit a high degree of pedogenic reorganization in the upper 0.7–0.8 m.

In the site of study, the soil is classified as Petrocalcic Argiudol (Soil Survey Staff 2014), and horizons sequence in afforested plot is (Figure 2(a)): (1) O (0–0.02 m), mostly composed of organic matter of litter in decomposition; (2) Ap (0.02–0.27 m), a clay loamy textured layer, with moderate and fine granular structure, very friable; (3) BA (0.27–0.47 m), clay loamy textured, with moderate subangular blocks that break to granular, very friable; (4) Bt1 (0.47–0.59 m), clay loamy textured with medium, moderate regular composite prisms that break into medium, weak and thin subangular prisms and abundant clay films; (5) Bt2 (0.59–0.76 m), with clay loamy texture, moderate, regular composite prisms that break into angular and sub-angular blocks, very hard and common clay-skin films, abrupt lower limit; and (6) 2CKkm (>0.76 m), a continuous and massive horizon nearly fully cemented by carbonates, very hard and practically impossible to penetrate by means of hand tools (petrocalcic horizon).

This study was developed in a 7.8 ha Eucalyptus viminalis forestation (Figure 1(b)), planted in 1999 with trees spacing of 2.5 by 3 m. However, 33% of trees are now dead due to the high density of planting and absence of management, such as thinning, fertilization and irrigation (remains ∼893 stems ha⁻¹). The undergrowth is practically absent, mostly Cynodon dactylon (except in the clearings, where it is less dominant), and species of asteraceous as Carduus acanthoides and Dispsacus fullonum to a lesser extent.

Grassland plot is adjacent to a forest one, and it covers an area of 3.6 ha (Figures 1(b) and 1(c)). In this plot, the dominant plant species in the form of mono specific tussocks were poaceous (Paspalum quadrifarium, Stipa brachychaeta, Pithochaetium spp., Mellica spp. and Bothriochloa laguroides) and apiaceous (Eryngium paniculatum). In addition, specimens of asteraceous, such as Carduus acanthoides, Dispsacus fullonum, Bellis spp., and Hypochaeris spp. were also identified.

Atmospheric parameters (temperature, relative humidity (RH), rainfall, etc.) were measured in both grassland and afforested plots. Thus, air temperature and RH were hourly recorded in the afforested plot with temperature and humidity Cavadevices sensors, and logger (models TC1047A, HIH4000, and 2K14, respectively; Cavadevices.com, Inc., CABA, Bs.As., AR). In the same plot, rainfall and stem flow were collected with a pluviometer and funnel installed on the tree trunk, respectively (monthly frequency).

As regards the atmospheric parameters in the grassland plot, only rainfall was monthly collected. Remaining parameters (i.e., temperature, RH and wind velocity) were obtained with an hourly frequency from two regional weather stations (see ‘La Chiquita’ and ‘La Germania’ stations in Figure 1(a)), located 4.3 and 8.2 km from the site of the study, respectively (BDH 2016).

Water table depth, electrical conductivity and temperature were hourly recorded using a CTD-10 Decagon sensor (Decagon Devices, Inc., Pullman, WA, USA) installed at 6 m depth in two boreholes in both forest and grassland plots, respectively (Figure 2(b)). These sensors use a vented differential pressure transducer to measure the pressure from the water column to determine the water depth. These measurements were monthly verified with a manual water level meter, whereas was verified in the laboratory from the collected groundwater samples.

To achieve a better exploration of the vertical soil profile, on 11/17/2016, two-dimensional electrical resistivity tomographies (ERT) cross sections in both grassland and afforested plots were collected (AGI Super String R1/IP equipment with 56 electrodes). Thus, electrodes were installed with a separation of 1 m, whereby electric current was injected and measured alternately and automatically to generate a grid of points of apparent resistivity. These values were inversed with the Earth Imager 2D code (Advanced Geosciences, Inc. 2005) to obtain the real resistivities (Dietrich et al. 2014), and then it was conversed to bulk electrical conductivities (⁠⁠, dS m⁻¹).

Water flow in the continuum phreatic aquifer–soil–atmosphere was simulated using the Min3P code (Mayer et al. 2002; Bea et al. 2012, 2016; Gamazo et al. 2015, see complete details of the numerical formulation in the Supplementary Material, SM. It solves the water flow under unsaturated conditions through the expanded Richards equation (Equation SM1; Saito et al. 2006; Bea et al. 2012), including vapor transport. The flow of liquid water is mainly a function of the physical properties of the soil, and the hydraulic and thermal gradients expressed through Darcy's law (Equation SM2, see details in Bea et al. 2012). To correctly solve the water flows in the soil, vadose zone implies to know the equilibrium between the liquid water and the gaseous phases, the latest one strongly temperature dependent. Thus, energy balance implemented in Min3P (Philip & De Vries 1957; Voss & Provost 2008; Bea et al. 2012) was applied in our analysis (Equation SM3). Both water flow and energy balances depend on climatic parameters, such as rainfall, RH and solar radiation. Physical evaporation, rainfall, runoff and forest water uptake are considered through a source/sink term in the flow equation (Equation SM5). The actual forest transpiration rate in Equation SM5 depends on soil moisture (Ghazanfari et al. 2015), and it is constrained by crop evapotranspiration under standard conditions ⁠. The last one is computed according to the Penman-Monteith formulation updated by Allen et al. (2006) and includes the crop coefficients (Kc) to E. viminalis extracted from Besteiro (2014). Thus, is computed as the sum of the water uptake contributions for each discrete soil volume representing the entire rhizosphere (Equations SM9–11).

To clarify further, monitoring results for grassland and forest plots are shown for comparison to assess the afforestation impact on water fluxes. Flow and transport modeling were only applied on the forest plot to keep the focus on the paper (i.e., methodology to estimate phreatophytes water consumption). Model domain for the afforested plot consisted of a 5.4 m soil column (Figure 2(b)), and it included the six identified soil horizons. Note that the organic matter horizon was excluded in the model due to its thickness (2 cm) since it was expected to have a low impact on water fluxes. Physical soil model parameters were obtained from sampling each horizon (see Table SM1 in SM). Soil particle size was obtained using the hydrometer method (Bouyoucos 1936), whereas the dry bulk density was measured in 169.65 cm³ undisturbed soil cores, after drying them for 24 h at 105 °C. Parameters for the van Genuchten model were estimated using the Rosetta Lite v.1.1 software (Schaap et al. 2001; Arrey et al. 2017) based on texture and bulk density (Table SM1 in SM). On the other hand, root length densities for each soil horizon were automatically calibrated coupling the Min3P code with PEST (Méndez et al. 2013; Doherty 2016), using the VWC (Ghazanfari et al. 2015) as an observation target. Calibration was constrained by the sap flow measurements in the time intervals when they were available.

The model is simulated in a 1D domain (i.e., horizontal water and heat exchanges were assumed to be zero), and recorded atmospheric boundary conditions were imposed at the top (Figure 2(b)). Recorded atmospheric temperature, RH, solar radiation and rainfall used in the simulations are shown in Figures 3(a) and 4(a). Note that hourly RH and rainfall records were obtained by combining data collected from two weather stations (‘La Chiquita’ and ‘La Germania’). They were combined by selecting one that best fit the in situ (monthly) precipitation records in each time interval. Thus, the atmospheric RH varied between 0.11 and 1, whereas the temperature varied between 38.1 and −5.0 °C from summer to winter, respectively. On the other hand, rainfall events increased their frequency in summer (rainy season), while they were scarce in winter. Note that the forest tree rainfall interception was calculated from rainfall measurements in both plots and stem flow in the afforested plot. Besides, the solar radiation was computed using the daily heliophany data from the weather station scaled according to the forest soil cover (this last one was obtained from image analysis using the Image J software, US National Institute of Health, Bethesda, Maryland, USA, 1997–2012). Lastly, the wind speed of about 2.75 m s⁻¹ was reduced to 50% in the afforested plot with respect to the grassland one according to the findings of other authors for eucalyptus (e.g., Stigter et al. 2002; Souza et al. 2010). The transient prescribed head was imposed at the bottom of the model domain (Figure 2(b)) to mimic the observed phreatic level fluctuations during the simulated period of time (see Figure 3(d)).

Additionally, conservative transport (electrical conductivity tracer, Equation SM12 in SM) was simulated to assess salt accumulation in these systems. Thus, soil, water table, and rainfall electrical conductivity measurements were imposed as initial and boundary conditions.

Temporal evolution of the VWC at 0.2 and 0.5 m depth is shown in Figures 3(b) and 3(c), respectively. VWC increased their values in response to rainfall events in both grassland and afforested plots. However, the VWC evolution was different according to the plot. As regards the grassland plot at 0.2 m (Figure 3(b)), it progressively increased their values (i.e., water storage increased), particularly during the rainy period (P1, from 10 December 2015 to 19 March 2016, Figure 3(b)). VWC reached a maximum value of about 0.3 in P1, and it maintained a quasi-constant value during the consecutive period (P2, from 19 March 2016 to 10 September 2016, Figure 3(b)) characterized by scarce rainfall events, lower temperatures and atmospheric demand (i.e., low ⁠). In the consecutive rainy period (P3, Figure 3(b)), VWC showed a similar behavior to P1. It followed a similar behavior at 0.5 m depth reaching maximum values of around 0.4 during the period with scarce rainfall events (P2, Figure 3(c)).

As it was expected in the afforested plot, VWC quickly responded to each rainfall event at 0.2 m depth, but they immediately decreased to the soil wilting point (Figure 3(b)). On the contrary, VWC at 0.5 m depth were practically unperturbed by the rainfall events (Figure 3(c)).

As regards the soil temperature evolution at 0.2 and 0.5 m depth, during the rainy periods, the day/night differences at 0.2 m were almost 2.5 and 1 °C in the grassland and afforested plots, respectively (see P1 and P3 periods, Figure 4(b)). However, these differences decreased during the winter to 0.3 and almost zero at 0.5 m depth in both plots (see P2 in Figure 4(b)). The average temperature in summer was around 5 °C higher in the grassland than in the afforested plot at both 0.2 and 0.5 m depth, whereas they were similar during the winter season (see P2 period in Figures 4(b) and 4(c)).

Temporal evolution of in the grassland and afforested plots is shown in Figure 5. In both plots at 0.2 m depth, they varied between 0.45 and 2.76 dS m⁻¹ during the studied period of time (Figure 5) and the average values in both plots were practically identical (1.55 and 1.82 dS m⁻¹ for the grassland and afforested plots, respectively). However, was slightly higher at 0.5 m depth in the grassland plot (range values between 0.14 and 6.45 dSm⁻¹) than in the afforested one (range values between 1.20 and 4.28 dS m⁻¹). Lastly, the average values at 0.5 m depth were 3.69 and 2.92 dS m⁻¹, in both plots, respectively (Figure 5).

The water table depth evolved in a similar way in both plots (Figure 3(d)), it started in the highest level at the beginning of the studied period of time (3.8 and 3.9 m of depth for the grassland and afforested plots, respectively), and it drew down almost 1 m during the first 125 days, in spite of the fact it corresponded to the first rainy period (see P1 period in Figure 3(d)). Then, it reached a slightly steady state during the P2 (Figure 3(d)) between 4.6 and 4.8 m and finally, it continued falling to 5.4 m depth. Water table fluctuations associated with tree transpiration were not observed in the afforested plot as it was noted in other sites of the Pampean region (e.g., see Engel et al. 2005; Jobbágy & Jackson 2007).

As regards the evolution of the phreatic aquifer temperature, it was stable along the period of time studied with negligible variations around 0.5 and 1 °C in both grassland and afforested plots, respectively (Figure 4). Concerning the evolution of the phreatic aquifer electric conductivities ⁠, their values were consistent with those ones measured for the regional groundwater system (∼0.66 dS m⁻¹; BDH 2016; Figure 5).

Through ERT, we obtained the two-dimensional cross sections in both grassland and afforested plots, respectively (Figure 6). Both and (soil 5TE Decagon sensors) varied according to a set of soil characteristics, such as moisture content (VWC), salinity, temperature and particle size.

In this case, assuming that mainly responded to soil VWC, it was possible to identify the water table position at 4.5 m depth (which coincided with CTD-10 sensor measurements, see Figure 3). On the other hand, petrocalcic horizon could be identified along both cross sections (slightly shallower in the grassland plot), but the forest cross one appeared to be more resistive (i.e., drier) than the grassland one above and below the petrocalcic horizon.

Forest transpiration rates computed from sap flow measurements (Figure 7(c)) were discontinued along the studied period of time due to an equipment malfunction. Only two-time intervals were assessed: (1) from 2 February 2016 to 30 March 2016 and (2) from 23 December 2016 to 5 January 2017. On the first-time interval, transpiration rates fluctuated between 1.95 and 4.87 mm day⁻¹ (daily average) with an average value of about 3.38 mm day⁻¹. On the second-time interval, measured values were 0.97, 5.38, and 2.69 mm day⁻¹ for minimum, maximum and average transpiration rates, respectively.

A base case and three alternative scenarios were simulated to assess a sensitivity analysis of the model parameters and potential conceptual models. Model parameters used in these simulations are detailed in Table SM1 (SM).

Simulated VWC and T showed, in general, a successful fitting with the measurements, as indicated by the global correlation indexes (see Table SM2 in SM). Simulated VWC at different depths increased as a consequence of the rainfall events and the quick forest water uptake until the soil wilting point (see results at 0.2 and 0.5 m depths in Figures 7(a) and 7(b), respectively). However, some differences between the measured and simulated results (see points 1 and 3 in Figure 7(a)) were predicted. On the other hand, note that the VWC are practically unperturbed at 0.5 m depth (Figure 7(b)). Modeling results also suggested that forest transpiration was an active process during both rainy periods covered in this study, always constrained by the atmospheric demand through (see points 1 and 3 in Figure 7(c)). The predicted mean transpiration rate was 2.08 mm day⁻¹ in the entire period of time (1.97 mm day⁻¹ in 2016, Table SM3), reaching values in the time intervals with sufficient available water stored (see point 1 in Figure 7(c)). Minimum transpiration rates of about 0.51 mm day⁻¹ (e.g., see point 2 in Figure 7(c)) were predicted when the VWC reached the soil wilting point (e.g., see point 2 in Figure 7(a)). Total simulated transpiration for the entire period of time was calculated with ∼918 mm (∼723 mm in 2016), with a ∼13% groundwater usage (i.e., consistent with a forest hydraulically disconnected from the aquifer). Simulated monthly water balance showed that forest transpiration rates practically depended on rainfalls (Figure 7(d)). In addition, computed transpiration rates were also consistent with the other ones estimated in the same region for E. camaldulensis (e.g., 250–500 mm yr⁻¹, 500–700 stems ha⁻¹; Engel et al. 2005) and E. viminalis (e.g., 922–976 mm yr⁻¹, 822–1,440 stems ha⁻¹; Besteiro 2014). Total vapor flow through the model domain was represented by E (evaporation) and condensation (Table SM3), which left a net water efflux of 71.46 mm yr⁻¹.

Simulated temperature evolution for the afforested plot at 0.2 and 0.5 m depth is shown in Figure SM2. As it was expected, they mainly responded to the incident solar radiation on the soil surface and to air temperature. Note that at 0.2 m depth temperature was approximately 2–3 °C higher than those ones measured during the rainy periods (P1 and P3), and these differences were minimal in the period of scarce rainfall events (P2). On the other hand, modeling results suggested a similar behavior for temperatures at 0.5 m depth (Figure SM2).

Three alternative scenarios were simulated to assess the model sensitivity of the chosen main parameters: (1) the root length density distribution (SC_RLD_UPPER, and SC_RLD_PHREATIC scenarios) and (2) the empirical parameters controlling the forest water uptake (SC_STRESSED, i.e., the crop coefficients and in Equation SM 10 (SM)). The results of these scenarios showed that the root length density calibrated by PEST was the most feasible scenario and the simulated VWC responded to this parameter even below the petrocalcic horizon (Figure SM3). In addition, the and crop coefficients were adequate over most of the simulated period of time, but some physiologic-stress responses not considered by the Min3P code caused some lesser important differences on simulated VWC in short periods of time (Figure SM4). More details of alternative scenarios are available in SM.

Water fluxes below afforested grassland with the presence of highly cemented horizons in the soil profile were perturbed as it was also observed in other sites (Jobbágy & Jackson 2004, 2007; Engel et al. 2005; Nosetto et al. 2008, 2009). However, these cemented horizons could imply a physical barrier to water flow and forest root growth, diminishing (and even preventing) their access to phreatic aquifer. In fact, in the site of study, water storage increased above this horizon in the grassland plot (mainly during the rainy periods), whereas it was immediately consumed by tree transpiration below the afforested plot, reaching the wilting point.

The lack of groundwater fluctuations, mainly attributed to forest transpiration in the site of study, hampered the assessment of forest groundwater consumption using a methodology based on them (eucalyptus, e.g., see Engel et al. (2005) and references therein). However, this fact could be explained in two ways, due to (1) a sort of hydraulic disconnection between the forest and the phreatic aquifer and (2) high aquifer hydraulic conductivities (∼5 m day⁻¹) that attenuated the perturbation induced by forest transpiration.

Process-based numerical modeling based on the modified Richards equation allowed us to assess here the forest transpiration rates and its component from groundwater (constrained by soil moisture measurements at different depths) in those systems with the aforementioned features. Modeling results suggested that the presence of petrocalcic horizons induced a partial hydraulic disconnection between the forest and the phreatic aquifer. They showed the occurrence of long water stress periods particularly in those time intervals of high atmospheric demand. During these stress periods, the forest was mainly feeding from rainwater stored above the petrocalcic horizon, with minimum access to phreatic aquifer. In fact, modeling results of a hypothetical forest hydraulically connected with the phreatic aquifer scenario overestimated the soil moisture above the petrocalcic horizon and did not predict realistic transpiration rates (close to ⁠, ∼1,200 mm yr⁻¹ in 2016).

Pore electrical conductivities ranged from 0.14 to 4.28 dS m⁻¹ in both grassland and afforested plots, respectively, and these values were a long way from affecting most of the crops’ productivity (e.g., see Maas 1985). Based on numerical results presented here, the potential risk of salinization of these soils with petrocalcic horizons was unlikely in the short and medium term (years or decades). Forest available water was mainly restricted to shallow soil portions, where the recharge occurred from rainfall with slightly electrical conductivities. It was consistent with the transport modeling results simulated here, where salt accumulation was not predicted at short and medium term (Figure 5).

As regards the model limitations, as it was stated above, the water uptake approach (Battaglia & Sands 1997) implemented in Min3P was more rigorous from a porous media perspective, but it was less accurate from a tree physiology point of view. However, this model limitation was also detected in other water stress response functions commonly used in physical soil modeling (e.g., Feddes et al. 1974; Šimůnek & Hopmans 2009; Verma et al. 2014; Cai et al. 2017).

Productive soil salinization induced by forest groundwater consumption emerged as a negative consequence of grassland afforestation in flatland areas. Estimation of the forest groundwater consumption was crucial assessing the soil salinization potential. However, the presence of petrocalcic soil horizons, common features in flatland areas around the world, prevented its estimation using traditional methodologies based on water table fluctuations. Thus, process-based numerical modeling constrained by soil water content measurements appeared here as an alternative. It was successfully applied to estimate the forest groundwater consumption in monitored afforested grassland plots with these features in the sub-humid Pampean region.

In the site of study, the numerical methodology proposed here allowed us to quantify water flows and uptake in afforested grasslands. Modeling results suggested a maximum of forest groundwater total transpiration usage of ∼13%, whereas in other works the eucalyptus consumed ∼50–70% (e.g., see Thorburn & Walker 1994; Cramer et al. 1999; Engel et al. 2005). However, note that soil salinization was not observed in this system, in spite of the fact it was detected in other sites inside the same region with similar or less eucalyptus water consumption (calculated here in a range of 522–722 mm yr⁻¹ compared to ∼250–500 mm yr⁻¹; Jobbágy & Jackson 2004, 2007). In this case, petrocalcic horizons constituted a water flow barrier, and the forest uptook water stored above the petrocalcic horizon because it could not freely access the deeper water resources.

We gratefully acknowledge the members of the Instituto de Hidrología de Llanuras ‘Eduardo Jorge Usunoff’ (IHLLA) for their field and laboratory measurements support. This work was funded by Research Grants from FONCYT-MINCyT (PICT-2013-1223 and PICT-2015-0744) and INTA (PNFOR 1104073), Argentina.

The Supplementary Data for this paper is available online at http://dx.doi.org/10.2166/hydro.2019.093.