Abstract

Grasslands are extensively distributed in flatland areas around the world, such as the Pampas in South America. It is one of the most economically productive landscapes and, as in other regions, they are being replaced by forests at increasing rates. Soil salinization emerges as a negative consequence associated with water deficit and forest groundwater consumption (∼250–500 mm/yr, in this region). An assessment of forest groundwater consumption is crucial for risk evaluation of soil salinization on flatland environments. For this aim, numerical modeling based on physical/biological processes and atmospheric boundary conditions was successfully applied in monitored grassland and afforested plots. Modeling results suggested a partial hydraulic disconnection between forest and phreatic aquifer due to the presence of petrocalcic horizons. Forest transpiration estimates were approximately 13% of total groundwater usage. Forest water consumption was then restricted to that soil portions above the petrocalcic horizons. Estimated forest transpiration rates (∼723 mm/yr) were similar to and even exceeded those reported in salinized sites with similar features. However, the risk of salinization of these soils was unlikely, because forest transpiration was restricted to the upper soil portions filled with fresh rainwater. The petrocalcic horizon retained water and prevented both deep drainage and, indirectly, soil salinization.

INTRODUCTION

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−1 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−1) 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−1, (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.

METHODOLOGY

Site description

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.

Figure 1

(a) Global distribution of soil inorganic carbon within areas with the presence of calcic or petrocalcic horizons (USDA-NRCS 2019). (b) Location of the site of study, and (c and d) details of main features of the installed grassland and afforested plots. Tree density is ∼893 stems ha−1 in the afforested plot.

Figure 1

(a) Global distribution of soil inorganic carbon within areas with the presence of calcic or petrocalcic horizons (USDA-NRCS 2019). (b) Location of the site of study, and (c and d) details of main features of the installed grassland and afforested plots. Tree density is ∼893 stems ha−1 in the afforested plot.

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).

Figure 2

(a) Vertical distribution of the different horizons identified in the soil profile and (b) simulated model domain and boundary conditions.

Figure 2

(a) Vertical distribution of the different horizons identified in the soil profile and (b) simulated model domain and boundary conditions.

Field data

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−1). 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).

On the other hand, soil and phreatic aquifer parameters were hourly recorded from 10 December 2015 to 23 February 2017. Soil volumetric water contents (VWC, m3 m−3), bulk electrical conductivity (, dS m−1) and temperature (T, °C) were measured using 5 TE Decagon sensors (Decagon Devices, Inc., Pullman, WA, USA) at two different depths above the petrocalcic horizon (0.2 and 0.5 m), by duplicate in both grassland and forest plots (Figure 2(b); Pumo et al. 2015). VWC measurements were validated comparing those ones obtained with a gravimetric method at the same depths (on 22 June 2016 and 21 December 2016 in the afforested plot and only on 22 June 2016 for the grassland plot). Based on the sensor measurements (which depend on VWC, solute concentrations and the soil particles conductivities, respectively), the pore electrical conductivity (, which depends mainly on solute concentrations in soil pore water) was computed according to Hilhorst (2000): 
formula
(1)
where is the permittivity of pore water, is the permittivity of bulk soil and is the permittivity of soil when bulk electrical conductivity is negligible. is calculated as a function of the soil temperature (, measured in the same probe) according to 
formula
(2)
where is measured by the sensor and is an empirical factor fixed here with a value of 4.1, as it was suggested by Hilhorst (2000).

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−1).

Sap flow corresponding to 9 trees was measured using Granier-type heat dissipation sensors (e.g., see Granier 1985, 1987; Figure 1(d)). Thus, two Granier probes were installed along each tree trunk at 1 m high, separated 0.1 m between them, using a socket punch to remove bark and phloem and a low speed drill to perforate the xylem. Only the upper probe was constantly heated, and the temperature difference between probes was measured and recorded every 30 minutes with an automatic data logger (Campbell CR10X, Campbell Scientific, Logan, Utah, USA). Thus, the sap velocity (, mm s−1) was calculated with the relationship stated by Granier (1985, 1987): 
formula
(3)
where k is related to the temperature difference between the two probes : 
formula
(4)
where is the temperature difference when .
Forest transpiration calculation implies to know the active xylem area of each measured tree. To determine the sapwood depth, wooden dowels were taken from several trees using a Pressleer borer (generally used by foresters). The sap flow rate (, mm s−1) was calculated according to 
formula
(5)
where is the sapwood and tree cross-sectional areas ratio (–). Temperature differentials were corrected for the proportion of heated probe in contact with sapwood following Clearwater et al. (1999). Finally, forest transpiration (, mm s−1) was calculated considering the basal area (, m2 ha−1) obtained with the Bitterlich method (Bitterlich 1948) as follows: 
formula
(6)

Modeling

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 cm3 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−1 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)).

Figure 3

Temporal evolution of (a) rainfall and RH, (b and c) soil water VWC in the grassland and afforested plots at 0.2 and 0.5 m depth, respectively, and (d) accumulated rainfall and water table depth (water table depth was not recorded in the grassland plot after September 2016 due to equipment malfunction).

Figure 3

Temporal evolution of (a) rainfall and RH, (b and c) soil water VWC in the grassland and afforested plots at 0.2 and 0.5 m depth, respectively, and (d) accumulated rainfall and water table depth (water table depth was not recorded in the grassland plot after September 2016 due to equipment malfunction).

Figure 4

Temporal evolution of (a) air temperatures and solar radiation, soil temperatures at (b) 0.2 and (c) 0.5 m depths, and (d) phreatic aquifer temperature below the grassland and afforested plots.

Figure 4

Temporal evolution of (a) air temperatures and solar radiation, soil temperatures at (b) 0.2 and (c) 0.5 m depths, and (d) phreatic aquifer temperature below the grassland and afforested plots.

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.

RESULTS

Field data

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−1 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−1 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−1) than in the afforested one (range values between 1.20 and 4.28 dS m−1). Lastly, the average values at 0.5 m depth were 3.69 and 2.92 dS m−1, in both plots, respectively (Figure 5).

Figure 5

Temporal evolution of soil pore electrical conductivity at (a) 0.2 and (b) 0.5 m depth (simulated electrical conductivity is also shown), and (c) the one corresponding to phreatic aquifer below the grassland and afforested plots.

Figure 5

Temporal evolution of soil pore electrical conductivity at (a) 0.2 and (b) 0.5 m depth (simulated electrical conductivity is also shown), and (c) the one corresponding to phreatic aquifer below the grassland and afforested plots.

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−1; 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.

Figure 6

ETR. Two-dimensional electrical conductivity transects across grassland and afforested plots, respectively, in the site of study. They were collected on 17 November 2016.

Figure 6

ETR. Two-dimensional electrical conductivity transects across grassland and afforested plots, respectively, in the site of study. They were collected on 17 November 2016.

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−1 (daily average) with an average value of about 3.38 mm day−1. On the second-time interval, measured values were 0.97, 5.38, and 2.69 mm day−1 for minimum, maximum and average transpiration rates, respectively.

Figure 7

Modeling results for the base case scenario. Temporal evolution for the simulated and measured VWC at (a) 0.2 m, (b) 0.5 m depth, (c) measured and simulated actual forest transpiration rates, and (d) monthly water balance.

Figure 7

Modeling results for the base case scenario. Temporal evolution for the simulated and measured VWC at (a) 0.2 m, (b) 0.5 m depth, (c) measured and simulated actual forest transpiration rates, and (d) monthly water balance.

Flow modeling results

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).

Base case simulation

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−1 in the entire period of time (1.97 mm day−1 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−1 (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−1, 500–700 stems ha−1; Engel et al. 2005) and E. viminalis (e.g., 922–976 mm yr−1, 822–1,440 stems ha−1; 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−1.

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).

Alternative simulated scenarios (model sensitivity)

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.

DISCUSSION

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−1) 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−1 in 2016).

Pore electrical conductivities ranged from 0.14 to 4.28 dS m−1 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).

CONCLUSIONS

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−1 compared to ∼250–500 mm yr−1; 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.

ACKNOWLEDGEMENT

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.

REFERENCES

REFERENCES
Advanced Geosciences
2005
EarthImager 2D, resistivity and IP inversion software, version 2.2.8. Instruction manual. Advanced Geosciences, Austin, TX
, p.
139
.
Allen
G. R.
,
Pereira
L. S.
,
Raes
D.
&
Smith
M.
2006
Crop Evapotranspiration- Guidelines for Computing Crop Water Requirements-FAO
.
Irrigation and drainage. Water Resources, Development and Management Service
,
Rome
,
Italy
.
1998, paper 56
, p.
323
.
Arrey
I. A.
,
Odiyo
J. O.
,
Makungo
R.
&
Kataka
M. O.
2017
Effect of hysteresis on water flow in the vadose zone under natural boundary conditions, Siloam Village case study, South Africa
.
Journal of Hydroinformatics
20
(
1
),
88
99
.
doi:10.2166/hydro.2017.091
.
Battaglia
M.
&
Sands
P.
1997
Modelling site productivity of Eucalyptus globulus in response to climatic and site factors
.
Functional Plant Biology
24
(
6
),
831
850
.
https://doi.org/10.1071/PP97065
.
BDH
2016
Base de Datos Hidrológicos
.
Available from: http://www.azul.bdh.org.ar (accessed March 2017)
.
Bea
S. A.
,
Wilson
S. A.
,
Mayer
K. U.
,
Dipple
G. M.
,
Power
I. M.
&
Gamazo
P.
2012
Reactive transport modeling of natural carbon sequestration in ultramafic mine tailings
.
Vadose Zone Journal
11
(
2
).
https://doi.org/10.2136/vzj2011.0053
.
Besteiro
S.
2014
Evaluación de la influencia hidrológica de forestaciones en la llanura pampeana. Tesis de Doctorado
,
Universidad Nacional de La Plata
,
La Plata, Buenos Aires
. .
Bitterlich
W.
1948
Die Winkelzahl probe (The angle count sample)
.
Allgemeine Forst- und Holzwirtschaftliche Zeitung
59
,
4
5
.
Bouyoucos
G. J.
1936
Directions for making mechanical analyses of soils by the hydrometer method
.
Michigan Agricultural & Experimental Station
269
,
225
229
.
Burgos
J. J.
1968
El clima de la Provincia de Buenos Aires en relación con la vegetación y el suelo
. In:
Flora de la Provincia de Buenos Aires
(
Cabrera
A. L.
ed.).
Parte 1
.
Colección Científica INTA
,
Buenos Aires
, pp.
33
99
.
Burgos
J. J.
&
Vidal
A. L.
1951
Los climas de la República Argentina según la nueva clasificación de Thornthwaite
.
Meteoros
1
,
1
32
.
Cai
G.
,
Vanderborght
J.
,
Couvreur
V.
,
Mboh
C. M.
&
Vereecken
H.
2017
Parameterization of root water uptake models considering dynamic root distributions and water uptake compensation
.
Vadose Zone Journal
.
https://doi.org/10.2136/vzj2016.12.0125
.
Calder
I. R.
1998
Water use by forests, limits and controls
.
Tree Physiology
18
(
8–9
),
625
631
.
https://doi.org/10.1093/treephys/18.8-9.625
.
Carreño
L.
,
Frank
F. C.
&
Viglizzo
E. F.
2012
Tradeoffs between economic and ecosystem services in Argentina during 50 years of land-use change
.
Agriculture, Ecosystems and Environment
154
,
68
77
.
https://doi.org/10.1016/j.agee.2011.05.019
.
Clearwater
M. J.
,
Meinzer
F. C.
,
Andrade
J. L.
,
Goldstein
G.
&
Holbrook
N. M.
1999
Potential errors in measurement of non uniform sap flow using heat dissipation probes
.
Tree Physiology
19
,
681
687
.
https://doi.org/10.1093/treephys/19.10.681
.
Cramer
V. A.
,
Thorburn
P. J.
&
Fraser
G. W.
1999
Transpiration and groundwater uptake from farm forest plots of Casuarina glauca and Eucalyptus camaldulensis in saline areas of southeast Queensland, Australia
.
Agricultural Water Management
39
(
2–3
),
187
204
.
https://doi.org/10.1016/S0378-3774(98)00078-X
.
Devito
K.
,
Mendoza
C.
&
Qualizza
C.
2012
Conceptualizing Water Movement in the Boreal Plains. Implications for Watershed Reconstruction
.
https://doi.org/10.1017/CBO9781107415324.004
.
Dietrich
S.
,
Weinzettel
P. A.
&
Varni
M.
2014
Infiltration and drainage analysis in a heterogeneous soil by electrical resistivity tomography
.
Soil Science Society of America Journal
78
(
4
),
1153
1167
.
https://doi.org/10.2136/sssaj2014.02.0062
.
Doherty
J. E.
2016
Model-independent Parameter Estimation User Manual Part I: Pest Sensan and Global Optimisers
.
Watermark Numerical Computing
,
Brisbane
,
Australia
, p.
390
.
Duniway
M. C.
,
Herrick
J. E.
&
Monger
H. C.
2010
Spatial and temporal variability of plant-available water in calcium carbonate-cemented soils and consequences for arid ecosystem resilience
.
Oecologia
163
(
1
),
215
226
.
https://doi.org/10.1007/s00442-009-1530-7
.
Engel
V.
,
Jobbágy
E. G.
,
Stieglitz
M.
,
Williams
M.
&
Jackson
R. B.
2005
Hydrological consequences of Eucalyptus afforestation in the Argentine Pampas
.
Water Resources Research
41
(
10
),
1
14
.
https://doi.org/10.1029/2004WR003761
.
Feddes
R. A.
,
Bresler
E.
&
Neuman
S. P.
1974
Field test of a modified numerical model for water uptake by root systems
.
Water Resources Research
10
(
6
),
1199
1206
.
https://doi.org/10.1029/WR010i006p01199
.
Gamazo
P.
,
Slooten
L. J.
,
Carrera
J.
,
Saaltink
M. W.
,
Bea
S.
&
Soler
J.
2015
PROOST: object-oriented approach to multiphase reactive transport modeling in porous media
.
Journal of Hydroinformatics
18
(
2
),
310
328
.
doi:10.2166/hydro.2015.126
.
Geary
T. F.
2001
Afforestation in Uruguay: Study of a changing landscape
.
Journal of Forestry
99
(
7
),
35
39
.
https://doi.org/10.1093/jof/99.7.35
.
Ghazanfari
S.
,
Pande
S.
,
Cheema
M. J. M.
,
Alizadeh
A.
&
Farid
A.
2015
The role of soil moisture accounting in estimation of soil evaporation and transpiration
.
Journal of Hydroinformatics
18
(
2
),
329
344
.
doi:10.2166/hydro.2015.114
.
Granier
A.
1987
Evaluation of transpiration in a Douglas-fir stand by means of sap flow measurements
.
Tree Physiology
3
(
4
),
309
320
.
https://doi.org/10.1093/treephys/3.4.309
.
Gyenge
J. E.
,
Fernández
M. E.
,
Salda
G. D.
&
Schlichter
T. M.
2002
Silvopastoral systems in Northwestern Patagonia II: water balance and water potential in a stand of Pinus ponderosa and native grassland
.
Agroforestry Systems
55
(
1
),
47
55
.
https://doi.org/10.1023/A:1020238330817
.
Herron
N.
,
Davis
R.
&
Jones
R.
2002
The effects of large-scale afforestation and climate change on water allocation in the Macquarie River catchment, NSW, Australia
.
Journal of Environmental Management
65
(
4
),
369
381
.
https://doi.org/10.1016/S0301-4797(02)90562-1
.
Heuperman
A.
1999
Hydraulic gradient reversal by trees in shallow water table areas and repercussions for the sustainability of tree-growing systems
.
Agricultural Water Management
39
(
2–3
),
153
167
.
https://doi.org/10.1016/S0378-3774(98)00076-6
.
Hilhorst
M. A.
2000
A pore water conductivity sensor
.
Soil Science Society of America Journal
64
(
6
),
1922
1925
.
https://doi.org/10.2136/sssaj2000.6461922x
.
Jackson
R. B.
,
Carpenter
S. R.
,
Dahm
C. N.
,
McKnight
D. M.
,
Naiman
R. J.
,
Postel
S. L.
&
Running
S. W.
2001
Water in a changing world
.
Ecological Applications
11
(
4
),
1027
1045
.
https://doi.org/10.1890/1051-0761(2001)011[1027:WIACW]2.0.CO;2
.
Jackson
R. B.
,
Banner
J. L.
,
Jobbágy
E. G.
,
Pockman
W. T.
&
Wall
D. H.
2002
Ecosystem carbon loss with woody plant invasion of grasslands
.
Nature
418
(
6898
),
623
.
http://dx.doi.org/10.1038/nature00910
.
Jobbágy
E. G.
&
Jackson
R. B.
2003
Patterns and mechanisms of soil acidification in the conversion of grasslands to forests
.
Biogeochemistry
64
(
2
),
205
229
.
https://doi.org/10.1023/A:1024985629259
.
Jobbágy
E. G.
&
Jackson
R. B.
2004
Groundwater use and salinization with grassland afforestation
.
Global Change Biology
10
(
8
),
1299
1312
.
https://doi.org/10.1111/j.1365-2486.2004.00806.x
.
Jobbágy
E. G.
&
Jackson
R. B.
2007
Groundwater and soil chemical changes under phreatophytic tree plantations
.
Journal of Geophysical Research: Biogeosciences
112
(
2
),
1
15
.
https://doi.org/10.1029/2006JG000246
.
Kuznetsova
A.
&
Khokhlova
O.
2015
Dynamics and genesis of calcic accumulations in soils and sediments of the Argentinean Pampa
.
International Journal of Sediment Research
30
(
3
),
179
189
.
https://doi.org/10.1016/j.ijsrc.2014.11.002
.
Maas
E. V.
1985
Crop tolerance to saline sprinkling water
.
Plant and Soil
89
(
1–3
),
273
284
.
https://doi.org/10.1007/BF02182247
.
Mayer
K. U.
,
Frind
E. O.
&
Blowes
D. W.
2002
Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions
.
Water Resources Research
38
(
9
),
1174
.
https://doi.org/10.1029/2001WR000862
.
Méndez
M.
,
Araya
J. A.
&
Sánchez
L. D.
2013
Automated parameter optimization of a water distribution system
.
Journal of Hydroinformatics
15
(
1
),
71
85
.
doi:10.2166/hydro.2012.028
.
Nosetto
M. D.
,
Jobbágy
E. G.
,
Tóth
T.
&
Di Bella
C. M.
2007
The effects of tree establishment on water and salt dynamics in naturally salt-affected grasslands
.
Oecologia
152
(
4
),
695
705
.
https://doi.org/10.1007/s00442-007-0694-2
.
Nosetto
M. D.
,
Jobbágy
E. G.
,
Tóth
T.
&
Jackson
R. B.
2008
Regional patterns and controls of ecosystem salinization with grassland afforestation along a rainfall gradient
.
Global Biogeochemical Cycles
22
(
2
),
1
12
.
https://doi.org/10.1029/2007GB003000
.
Nosetto
M. D.
,
Jobbágy
E. G.
,
Jackson
R. B.
&
Sznaider
G. A.
2009
Reciprocal influence of crops and shallow ground water in sandy landscapes of the Inland Pampas
.
Field Crops Research
113
(
2
),
138
148
.
https://doi.org/10.1016/j.fcr.2009.04.016
.
Nosetto
M. D.
,
Jobbágy
E. G.
,
Brizuela
A. B.
&
Jackson
R. B.
2012
The hydrologic consequences of land cover change in central Argentina
.
Agriculture, Ecosystems and Environment
154
,
2
11
.
https://doi.org/10.1016/j.agee.2011.01.008
.
Paruelo
J. M.
,
Guerschman
J. P.
,
Piñeiro
G.
,
Jobbagy
E. G.
,
Verón
S. R.
,
Baldi
G.
&
Baeza
S.
2006
Cambios en el uso de la tierra en Argentina y Uruguay: marcos conceptuales para su análisis
.
Agrociencia
10
(
2
),
47
61
.
https://doi.org/10.2307/2577037
.
Pazos
S. M.
&
Mestelan
A. S.
2002
Variability of depth to tosca in udoll sand soil classification, Buenos Aires Province, Argentina
.
Soil Science Society of America Journal
66
,
1256
1264
.
https://doi.org/10.2136/sssaj2002.1256
.
Philip
J.
&
De Vries
D.
1957
Moisture movement in porous materials under temperature gradients
.
Transactions, American Geophysical Union
38
(
2
),
222
232
.
https://doi.org/10.1029/TR038i002p00222
.
Prieto
A. R.
1996
Late quaternary vegetational and climatic changes in the Pampa Grassland of Argentina
.
Quaternary Research
45
(
1
),
73
88
.
https://doi.org/10.1006/qres.1996.0007
.
Pumo
D.
,
Francipane
A.
,
Lo Conti
F.
,
Arnone
E.
,
Bitonto
P.
,
Viola
F.
&
Noto
L. V.
2015
The SESAMO early warning system for rainfall-triggered landslides
.
Journal of Hydroinformatics
18
(
2
),
256
276
.
doi:10.2166/hydro.2015.060
.
Richardson
D. M.
1998
Forestry trees as invasive aliens
.
Conservation Biology
12
(
1
),
18
26
.
http://dx.doi.org/10.1111/j.1523-1739.1998.96392.x
.
Saito
H.
,
Simunek
J.
&
Mohanty
B.
2006
Numerical analysis of coupled water, vapor, and heat transport in the vadose zone
.
Vadose Zone Journal
5
,
784
800
.
https://doi.org/10.2136/vzj2006.0007
.
Salama
R. B.
,
Bartle
G. A.
&
Farrington
P.
1994
Water use of plantation Eucalyptus camaldulensis estimated by groundwater hydrograph separation techniques and heat pulse method
.
Journal of Hydrology
156
(
1–4
),
163
180
.
https://doi.org/10.1016/0022-1694(94)90076-0
.
Schaap
M. G.
,
Leij
F. J.
&
Van Genuchten
M. T.
2001
Rosetta: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions
.
Journal of Hydrology
251
(
3–4
),
163
176
.
https://doi.org/10.1016/S0022-1694(01)00466-8
.
Shreve
F.
&
Mallery
T. D.
1932
The relations of caliche to desert plants
.
Soil Science
35
,
99
113
.
Šimůnek
J.
&
Hopmans
J. W.
2009
Modeling compensated root water and nutrient uptake
.
Ecological Modelling
220
(
4
),
505
521
.
https://doi.org/10.1016/j.ecolmodel.2008.11.004
.
Soil Survey Staff
2014
Keys to Soil Taxonomy
(12th edn.)
United States Department of Agriculture. Natural Resources Conservation Service
,
Washington, DC
.
Souza
W. D.
,
Barbosa
O. R.
,
Marques
J. D. A.
,
Costa
M. A. T.
,
Gasparino
E.
&
Limberger
E.
2010
Microclimate in silvipastoral systems with eucalyptus in rank with different heights
.
Revista Brasileira de Zootecnia
39
(
3
),
685
694
.
https://doi.org/10.1590/S1516-35982010000300030
.
Stigter
C. J.
,
Mohammed
A. E.
,
Nasr Al-Amin
N. K.
,
Onyewotu
L. O. Z.
,
Oteng̀i
S. B. B.
&
Kainkwa
R. M. R.
2002
Agroforestry solutions to some African wind problems
.
Journal of Wind Engineering and Industrial Aerodynamics
90
(
10
),
1101
1114
.
https://doi.org/10.1016/S0167-6105(02)00224-6
.
Thorburn
P. J.
1997
Land management impacts on evaporation from shallow, saline water tables
. In:
Subsurface Hydrological Responses to Land Cover and Land Use Changes
. (T. Makoto, ed.),
Springer
,
Boston, MA, USA
. pp.
21
34
.
http://doi.org/10.1007/978-1-4615-6141-5_2
.
Thorburn
P. J.
&
Walker
G. R.
1994
Variations in stream water-uptake by Eucalyptus-camaldulensis with differing access to stream water
.
Oecologia
100
(
3
),
293
301
.
http://doi.org/10.1007/BF00316957
.
USDA-NRCS
2019
World Soil Resources
. Available from:
http://soils.usda.gov (accessed May 2019)
.
Van Genuchten
M. T.
1980
A closed-form equation for predicting the hydraulic conductivity of unsaturated soils
.
Soil Science Society of America Journal
44
(
5
),
892
898
.
http://doi.org/10.2136/sssaj1980.03615995004400050002x
.
Varni
M. R.
&
Usunoff
E. J.
1999
Simulation of regional-scale groundwater flow in the Azul River basin, Buenos Aires Province, Argentina
.
Hydrogeology Journal
7
(
2
),
180
187
.
https://doi.org/10.1007/s100400050190
.
Verma
P.
,
Loheide
S. P.
II.
,
Eamus
D.
&
Daly
E.
2014
Root water compensation sustains transpiration rates in an Australian woodland
.
Advances in Water Resources
74
,
91
101
.
http://dx.doi.org/10.1016/j.advwatres.2014.08.013
.
Voss
C.
&
Provost
A.
2008
SUTRA – A Model for Saturated-Unsaturated Variable-Density Ground-Water Flow with Solute or Energy Transport (Version 2.1)
.
USGS Water Resources Investigations Report, 2002–4231
.
Wright
J. A.
,
DiNicola
A.
&
Gaitan
E.
2000
Latin American forest plantations: opportunities for carbon sequestration, economic development, and financial returns
.
Journal of Forestry
98
(
9
),
20
23
.
https://doi.org/10.1093/jof/98.9.20
.
Zalba
S. M.
&
Villamil
C. B.
2002
Woody plant invasion in relictual grassland
.
Biological Invasions
4
(
1–2
),
55
72
.
https://doi.org/10.1016/S0167-8809(01)00209-2
.
Zárate
M. A.
2003
Loess of Southern South America
.
Quaternary Science Reviews
22
(
18–19
),
1987
2006
.
https://doi.org/10.1016/S0277-3791(03)00165-3
.
Zárate
M.
&
Mehl
A.
2010
Geología y geomorfología de la cuenca del arroyo del Azul, provincia de Buenos Aires, Argentina
. In
I Congreso Internacional de Hidrología de Llanura-Azul
,
21 Al 24 September 2010
,
Buenos Aires, Argentina
, pp.
81
91
.

Supplementary data