CC BY-NC-ND
0GEODERjVIA
Geoderma xxx (xxxx) xxx-xxx
Contents lists available at ScienceDirect
Geoderma
journal homepage: www.elsevier.com/locate/geoderma
Soil variability in La Violada Irrigation District (Spain): II Characterizing hydrologic and salinity features
M.T. Jiménez-Aguirrea'*, D. Isidoroa, A. Usónb
a Centro de Investigación y Tecnología Agroalimentaria de Aragón (CITA), Unidad de Suelos y Riegos (Unidad Asociada EEAD-CSIC), Avenida de Montañana 930, CP 50059 Zaragoza, Spain
b Universidad de Zaragoza, Escuela Politécnica Superior de Huesca, Departamento de Ciencias Agrarias y del Medio Natural, Carretera Cuarte s/n, CP 22071 Huesca, Spain
ABSTRACT
The recent modernization of 1.1 Mha of irrigated land in Spain calls for the evaluation of these transformations in terms of environmental impact and resource use efficiency. The available data for this evaluation has increased with the transformation (better, digital and spatially distributed data) allowing for the use of distributed soil water and solute movement models. But most hydrological models require soil hydrologic properties that are costly and time-consuming to gather and soil information in Spain is generally scarce. This paper focuses in analyzing the soil hydrologic features in La Violada Irrigation District (VID; a 5234 ha semi-arid irrigated area recently modernized in northeast Spain) usable in soil water models for the evaluation of the new irrigation system.
The recent soil map of the VID (presented in a companion paper) gathered the hydrologic and salinity properties of the horizons in the described soil units. The hydraulic conductivity (Ks) of the horizons was also assessed by the inverse auger-hole method. From these data, the VID was disaggregated in three homogeneous units according to their hydrologic features and Pedotransfer Functions (PTFs) were built for the whole VID (General Model) and separately for the homogeneous soil units (Distributed Model). These PTFs allowed for obtaining field capacity (FC) and wilting point (WP) from texture and organic matter, while Ks depended upon texture and gypsum content.
Apparently, there were no salinity issues in VID soils due to irrigation. The high Ca2+ and Mg2+ levels in the saturation extract resulted in generally low SAR, what along with the high gypsum and carbonate contents may help to prevent soil degradation by sodicity.
As a result, the homogeneous hydrologic zones defined in VID may be used to recommend specific irrigation practices and as the basis for the application of distributed soil water movement models. These hydrologic properties may be applied directly as inputs to the models while the PTFs may allow for setting adequate parameters in nearby areas with similar soils from more readily available soil information (texture, organic matter and gypsum).
ELSEVIER
ARTICLE INFO
Keywords:
Pedotransfer Functions Available Water Content Hydraulic conductivity Soil water balances Salinity
Gypsum rich soils
1. Introduction
The first part of this work presented the Violada Irrigation District (VID) soil map and the distribution of the main soil properties in depth and along the VID (thematic maps; Jimenez-Aguirre et al., 2017). Hydraulic conductivity, infiltration and soil water retention are the
essential input data for soil water movement models (Minasny and Hartemink, 2011; Nguyen et al., 2015); but the measurement of these properties is a generally complex and time-consuming process (Wagner et al., 2001; Wösten et al., 2001).
Pedotransfer Functions (PTF) provide an alternative to estimate hydrologic soil data [such as Field Capacity (FC), Wilting Point (WP) or
Abbreviations: ANOVA, Analysis of the Variance; AWC, Available Water Capacity; AWUA, Almudevar Water User Association; BD, Bulk Density; CCE, Calcium Carbonate Equivalent; CE, Coarse Elements; dfE, degrees of freedom of the Error; dfT, degrees of freedom of the Total; ECi:5, Electric Conductivity of the 1:5 soil-water extract; ECe, Electric Conductivity of the saturation extract; FC, Field Capacity; GC, Gypsum Content; IRF, Irrigation Return Flows; Ks, Saturated Hydraulic Conductivity; LSD, Least Significant Difference; n, number of submodels from the Distributed Model; N, number of samples; OM, Organic Matter; PTF, Pedotransfer Functions; RMSE, Root Mean Square Error; SAR, Sodium Absorption Ratio; SSE, Sum of Squares of the Error; SST, Sum Squares of the Total; VID, Violada Irrigation District; WP, Wilting Point; WUA, Water User Association * Corresponding author.
E-mail addresses: mtjimeneza@cita-aragon.es (M.T. Jimenez-Aguirre), disidoro@aragon.es (D. Isidoro), mauson@unizar.es (A. Uson). http://dx.doi.org/10.10167j.geoderma.2017.04.024
Received 21 October 2016; Received in revised form 7 April 2017; Accepted 26 April 2017
0016-7061/ © 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.Org/licenses/BY-NC-ND/4.0/).
Please cite this article as: Jiménez-Aguirre, M., Geoderma (2017), http://dx.doi.Org/10.1016/j.geoderma.2017.04.024
M.T. Jiménez-Aguirre et ai
Fig. 1. Map of La Violada Irrigation District: location, irrigation canals, soil map, and observation points [hydraulic conductivity tests (labelled "K"), soil pits ("C") and auger holes (x)] throughout the Violada Irrigation District [from Jimenez-Aguirre et al. (2017)].
Saturated Hydraulic Conductivity (Ks)] from more usual soil survey data as Texture, Organic Matter (OM) or Gypsum Content (GC) (Bouma, 1989; Wösten et al.; 1999, 2001). On the other hand, hydraulic properties, and therefore PTFs, present a temporal and spatial variability with a significant effect on the model results; thus calling for establishing specific, more accurate PTFs for areas with similar soil hydrologic characteristics (Distributed Model) (Franzmeier, 1991; Pachepsky et al., 2006; Wessolek et al., 2011; Wösten et al., 2001).
In Spain, the need for evaluation of the irrigation districts (around 1.1 Mha) recently modernized by the two Nation Irrigation Modernization Plans (MARM, 2002, 2006) is growing. A recently modernized district (2008-09), La Violada Irrigation District (VID), was selected to analyze the impacts of the conversion of a traditional surface irrigation system into sprinkler irrigation.
The VID (5234 ha) is located in the middle Ebro River Basin (Spain), a semi-arid region. VID has been widely studied since the 1980s as described in the companion paper, in regard to agricultural management (crops, irrigation and fertilization) and salts and nutrients loads in the irrigation return flows under traditional irrigation. The new irrigation system in VID, with detailed information about volume and schedule of water applied at hydrant level (Stambouli, 2012) has increased the available information, not only about irrigation, but about fertilization management and crop distribution as well. This calls for the application of distributed soil water models in VID that could yield better estimates of salts and nutrient removal from the district and allow for simulating the effects of plausible future scenarios. And in order to apply such models, better soil information is required.
Developing hydrologic PTFs (FC, WP, Ks) for the whole VID and for smaller homogeneous units with similar hydrological properties, may allow for defining homogeneous hydrological response units, the basis for applying distributed water balance models. These may improve the
studies performed in the VID about irrigation return flows, salt and pollutant loads, and water use. In addition, the use of Geographic Information Systems (GIS) permits a higher spatial disaggregation to apply all this information and run models.
Furthermore, the Ebro River Valley is highly vulnerable to saliniza-tion (Ibarra, 2004). Salinity is a persistent problem in many irrigated lands in arid and semi-arid regions (Díaz and Herrero, 1992). The risk of salinization or sodification should be checked in a semi-arid irrigated, gypsum-rich area as VID. The secondary salinization affecting irrigated areas has been described (Szabolcs and Várallyay, 1979) using the electrical conductivity of the saturated paste extract (ECe) or the 1:5 soil-water solution (EC1:5) indistinctly as salinity indicators [threshold ECe > 4 dS/m (United States Salinity Laboratory Staff -USSLS, 1954)] for salinity affected soils. However several authors reported the weak correlation between ECe and EC1:5 in soils where GC varies (Aragüés et al., 1986; Casby-Horton et al., 2015; Herrero and Bercero, 1991; Herrero et al., 2009, Moret-Fernández and Herrero, 2015; Nogués et al., 2006).
This work is the second part of an integrated study with the objective of characterizing soil variability (in regard to hydrologic and irrigation related properties) within the VID. Taking as a basis the soil units defined in a companion paper, the specific objectives of this paper are: (i) to define hydrologically homogenous soil zones in the VID; (ii) to define Pedotransfer Functions to link the soil map with soil hydrologic features (FC, WP and Ks) for the homogeneous zones; (iii) to analyze the salinity issues and the risk of salinization in the VID.
The results of this two-part study (combined with distributed irrigation data) could be applied in further modeling works requiring detailed input data or in decision-making processes at water user association (the VID or other irrigation districts with similar hydro-logical characteristics and scarce soil information) or higher level (such
as basin authority).
2. Site description
The Violada Irrigation District (VID; 5234 ha) is located in the middle Ebro River Basin in Spain (Fig. 1), in the medium reaches of La Violada Gully Basin upstream of the gauging station of La Pardina (no 230) of the Ebro River Basin Authority. The climate, geology and irrigation system of the VID are described in a companion paper (Jimenez-Aguirre et al., 2017), along with the previous soil, irrigation and environmental studies in the area.
2.1. The VID soil map
The first part of this work presented a semi-detailed soil map according to the Soil Taxonomy classification (Soil Survey Staff, 2014), focused on the genesis and variability of the VID soils linked to irrigation considerations (Jimenez-Aguirre et al., 2017). The resulting map had 13 soil units by the combination of five subgroups (Fig. 1): Typic Calcixerept (A), Petrocalcic Calcixerept (B), Gypsic Haploxerept (C), Typic Xerorthent (D) and Typic Xerofluvent (E); and six particle size families: Fine (1), Fine Silty (2), Fine Loamy (3), Coarse Loamy (4), Loamy (shallow soils) (5) and Loamy-skeletal (6).
The soil map was based on 34 soil pit descriptions (environmental and physiographic characteristics) and their laboratory determinations: texture, Calcium Carbonate Equivalent (CCE), GC, OM, FC, WP and bulk density (BD); and 33 auger-holes, as described in the companion paper.
3. Materials and methods 3.1. Soil data collection
The data from the soil map of the VID was employed and extended. The texture was determined by the pipette method (Soil Survey Staff, 2011) into clay (particles with diameter < 0.002 mm), silt (between 0.002 and 0.05 mm) and sand (0.05-2 mm) fractions. In this second part the silt fraction was discerned between fine silt (0.002-0.02 mm) and coarse silt (0.02-0.05 mm). The water content [determined at FC (- 33 kPa) and WP (- 1500 kPa) in % weight], CCE, GC and OM were determined as explained in the companion paper (Jimenez-Aguirre et al., 2017).
The pH was measured in a 1:2.5 soil:water solution with a pH meter GLP22 and the soil electrical conductivity in a 1:5 soil:water solution (EC1:5) with a Orion 5 Star conductivity meter. The saturated paste extract was prepared for the 110 horizon samples taken and analyzed for electrical conductivity (ECe), pH (pHe), main ions [cations (Na+, K+, Ca2+, Mg2+) and anions (Cl- and HCO3-)] and Sodium Absorption Ratio (SAR) following United States Salinity Laboratory Staff (USSLS, 1954). Cations were measured using a continuous flow analyzer (AutoAnalyzer3, Bran + Luebbe), Chloride with a chloride analyzer (Sherwood model 926), and bicarbonates by potentiometric titration. Sulfate ion (SO42 - ) was determined as the difference between cations and anions.
The Available Water Capacity (AWC) was determined as the difference between FC and WP (mm) down to 60 cm (or impervious horizon). This depth was chosen because 80% of the pit profiles presented "frequent" roots under the criteria of SINEDARES (CBDSA, 1983) only down to 60 cm. Deeper soil layers were not considered due to the scarce presence of roots. The AWC for A horizons and B horizons (down to 60 cm) was considered separately as well. The water contents at FC and WP (in mm) were calculated from the data in the companion paper: horizon thickness, BD and coarse elements (CE), along with FC and WP in % weight. All data came from the 34 soil pits and the 33 auger holes opened in the VID; Fig. 1; Jimenez-Aguirre et al., 2017). The CE in horizon layers was estimated visually as percent volume in
the pits.
From May to June 2013, Saturated Hydraulic Conductivity (Ks) was measured in 15 selected locations close to the soil pits considered most representative of each soil unit and labelled with the letter "K" (Fig. 1), using the inversed auger-hole method (Oosterbaan and Nijland, 1994). The auger holes were opened by manual drill down to a given depth according to the depth of the tested horizon in the representative pit. Two or three repetitions per horizon A and B (C if possible) were made. One of the 15 Ks tests made (the K-7 test) was deleted. This test was performed in a field employed for stacking work materials, crop residues and manure, possibly causing soil compaction. As a result the Ks values from K-7 were extremely low and considered unrealistic. The horizons were classified by their Ks following the Soil Survey Division Staff (1993).
3.2. Statistical procedures
The differences in FC, WP and AWC (FC minus WP) of the A and B horizons down to 60 cm between the subgroups were tested by means of an ANOVA based on the Least Significant Difference (LSD) test. The Ks was also tested for differences between the subgroups, horizons and particle size families. The limited number of observations did not allow for a multiple ANOVA; thus the three factors were taken as single factors. A t-test for paired samples was performed to compare the Ks from the tests with the values offered by the ROSSETA v1.1 software (Schaap et al., 2001) provided by the HYDRUS-1D v4 software (Simunek et al., 2008) using the texture and BD as input.
The Ks values of each horizon in each soil unit were assigned to the corresponding horizon in the soil map, so that the Ks could be linked to the other horizon data in the soil map by means of Pedotransfer Functions (PTF).
Pedotransfer Functions (PTF) for FC, WP (percent weight) and Ks (m/d) were tested upon the variables clay (or sand), silt (also fine silt), OM and GC in each horizon by means of a multiple regression, using all the observations available in VID (General Model). As Ks and GC showed a non-normal distribution (as assessed through the Kolmogorov-Smirnov test) the multiple regressions were performed on their natural logarithms (LnKs and LnGC) which followed a normal distribution. Only the significant variables were retained in the final regressions and the residuals of these final models were checked for normality (Kolmogorov-Smirnov test).
The residuals from the PTFs for FC and WP including the C horizon showed a non-normal distribution, thus only A and B horizons were employed for the regressions of FC and WP (with normal residuals). Altogether, 88 soil samples for FC and WP and 69 repetitions for Ks were used.
With the aim of disaggregating the General Model (obtaining PTFs for the soil subgroups or their combinations), the same PTFs were tested on 5 different combinations of the soil subgroups to define n homogeneous zones (2 or 3) in VID with similar hydrologic characteristics following different grouping criteria: "Calcixerepts" (subgroup combination: AB vs CDE), "Soil Order" (ABC-DE), "Geomorphology" (AB-CE-D), "Carbonates or Gypsum accumulations" (AB-C-DE) or "Particle size family" (123-456). The PTFs for FC and WP were calculated for these combinations and assessed by two statistical indicators: the adjusted coefficient of determination (R^dj) and the root mean square error (RMSE). In each regression, only the significant variables were retained and their residuals were also tested for normality (Kolmogorov-Smirnov). Higher R2adj and lower RMSE were the criteria to select the best fitting combination, the RMSE decrease being the main statistical indicator versus the R2adj increase. The best fitting combination was selected as the Distributed Model of the VID (sum of the submodels from the n homogeneous zones). The concordance correlation coefficient (Lin, 1989) between observed and estimated values was also calculated for each best fitting model to assess possible systematic errors.
Fig. 2. Two families of soil samples identified in the Violada Irrigation District by their gypsum content (GC).
The R2dj for the Distributed Model was calculated from the statistical parameters of the n submodels: sum of squares of the error (SSE), degrees of freedom of the error (dfE), sum squares of the total submodel (SST) and degrees of freedom of the total submodel (dfT); following the expression:
m=issEi
R2a4 = 1- -
11=1 f
11=1 sst
11=1 fi
Similarly, a General Model and a Distributed Model for the Ks were set up with PTFs by multiple regressions. For the Distributed Model the same three previous homogeneous zones were used for the fine textured soils. But for the coarser textures [Loamy (shallow) and Loamy-skeletal] a separate regression was used to establish a single PTF. After identifying two families of GC present in VID (Fig. 2) two kind of PTFs were tested: one with Quantitative GC and the other with Qualitative GC as variables. The former used LnGC as variable. The Qualitative GC was introduced as a dummy variable (value 1 for GC higher than 4%; value 0 for GC < 4%) which will only affect the constant of the multiple regression.
To evaluate the salinity in the VID, the EC1:5 was compared with the ECe for the 110 horizon samples. A principal components analysis was performed using the two standardized ECs as variables. Then, a cluster analysis (Ward's method on the standardized variables with the Euclidean distance) was performed on the EC principal components (EC-Cluster). Also, SAR was related to ECe and EC1:5 by means of a multiple regression analysis.
Finally a Factor Analysis was completed on the 110 horizon samples using as variables the CCE, GC, ECe; and cations Na + , Ca2 + and Mg2 + (K+ was excluded as it was lower than 5% in all samples in the VID; Fig. 3) and anions Cl- and HCO3- in the saturation extract. CCE and GC were selected for this analysis due to the importance of these variables
Fig. 3. Frequency of the main cations in the saturated paste extracts of the soil horizons sampled in the Violada Irrigation District expressed as percentage of the sum of the main cations Na+, K+, Ca2+ and Mg2+.
in describing the variability of the soil properties in the VID (Jimenez-Aguirre et al., 2017). The variables were standardized and the main principal components selected were rotated by the varimax method to achieve orthogonal factors (Harman, 1967).
4. Results and discussion
4.1. Available Water Capacity and Hydraulic Conductivity
The ANOVA analysis for FC, WP and AWC (in mm) for the A horizons showed differences between the subgroups (Table 1) with the AWC being higher for Gypsic Haploxerept and Typic Xerofluvent (C and E subgroups) than for the Typic and Petrocalcic Calcixerept and Typic Xerorthent, (A, B and D). B horizons showed slight differences close to significance (p < 0.05). Down to 60 cm (A and B horizons together) there were significant differences between the subgroups; with the Petrocalcic Calcixerept (B) having the lowest AWC followed by the Typic Calcixerept (A), and with the highest AWC found in the subgroups C, D and E, with no significant differences between them.
The range of AWC (in mm for the whole profile, 60 cm) obtained was similar to that proposed by Isidoro (1999) (down to 120 cm or impenetrable layer) for 5 soil classes based on the observations of Slatni (1996) although lower AWC (mm) values were expected due to the difference in the soil thickness considered. Slatni (1996) may have overestimated the water content drying the soil samples at 105 °C instead 60 °C as proposed later by Artieda et al. (2006) and Herrero et al. (2009) for gypsum-rich soils (to avoid counting gypsum constitution water as soil water). So the AWC found by Slatni (1996) in percentage of water content (AWCSlatni = 9.2%; difference between FC and WP in % weigh: average FC = 28.3% and WP = 19.1% weight) was lower than the AWC in this work (AWC = 11.4%; average FC = 24.8% and WP = 13.4% weight).
Although there was no difference between the thickness of the A and B horizon for the subgroups, the contribution of the B horizon to the total available water was higher in the Gypsic Haploxerept Typic Xerorthent and Typic Xerofluvent (subgroups C, D, and E; Valley bottoms) than in the Typic and Petrocalcic Calcixerept (subgroups A and B; Glacis).
The ANOVA test for Ks showed no significant differences between horizons (A and B) or subgroups; but there were significant differences in regard to texture and qualitative gypsum. The number of observations for Ks within each soil subgroup or particle size family did not allow for an ANOVA test on the differences between subgroups and families and other variables, but some clues might be inferred from the mean values (Table 2): for fine textured soils (1 to 4 particle size families), the Ks within a subgroup generally increased as the particle size increased, being highest for the Coarse Loamy family; while for the coarse textured soils [Loamy (shallow) and Loamy-skeletal particle size families], the Ks was somewhat lower, in the range of the Fine Silty and Fine Loamy particle size families (possibly suggesting that the Ks in the coarse particle size families is controlled by the fine fraction, loamy in classes 5 and 6; Table 2).
In relation to gypsum (taken as a qualitative variable), the lnKs was significantly higher for soils with gypsum (Ks = 0.39 m/d if GC > 4%) than for soils without gypsum (Ks = 0.19 m/d if GC > 4%), also when the coarse textured Loamy (shallow) and Loamy-skeletal families were excluded.
4.2. Pedotransfer Functions
PTFs were established for the Global Model of the VID and for the Distributed Model. Three homogenous zones (Fig. 4) with similar hydrologic characteristics resulted from the disaggregation process of the VID with the "Carbonates or Gypsum accumulations" criteria as the best fitted combination: Zone I (CaCO3 accumulations) corresponds with the Typic Calcixerept and Petrocalcic Calcixerept (A and B
Table 1
Average Field Capacity, Wilting Point and Available Water Capacity (mm) for the A horizons, B horizons (from A horizon limit to down to 60 cm) and whole profile down to 60 cm for the different subgroups in the VID. Different letters show significant differences between subgroups (within each horizon and to 60 cm deep) by the LSD test (P < 0.05).
Field Capacity (mm) Wilting Point (mm) Available Water Capacity (mm)
Subgroup A horizon B horizon 60 cm depth A horizon B horizon 60 cm depth A horizon B horizon 60 cm depth
A 104.5a 58.5a 163.0a 52.8a 24.3a 77.1a 51.6a 34.2ab 85.8b
B 110.8ab 16.7a 128.8a 60.7ab 8.8a 80.6ab 50.2a 8.0a 54.1a
C 167.8c 74.5ab 242.3b 94.3c 41.6ab 136.0c 73.5b 32.9a 106.3c
D 117.7a 101.7b 218.5b 60.8ab 53.1b 114.5bc 56.9a 48.6b 102.3c
E 156.6bc 81.0ab 243.3b 83.0bc 44.7ab 131.8c 73.6b 36.3ab 105.8c
subgroups). Zone II (Gypsum accumulations) comprised the Gypsic Haploxerept subgroup (C) and Zone III (superficial gypsum outcrops or very deep gypsum accumulations) included the Typic Xerorthent and Typic Xerofluvent subgroups (D and E). The definition of these zones were similar to nearby study areas (Nogues and Herrero 2003; Nogues et al., 2006). Separated PTFs were established in each of the three zones for the FC and WP. Distinct PTFs for the Ks were established for fine textured soils with the same three zones and a differentiated PTF for coarse textured zone with the Loamy (shallow) and Loamy-skeletal particle size families (5 and 6) (Fig. 4).
4.2.1. Field Capacity and Wilting Point
The Pedotransfer Functions in VID linked FC and WP (% weight) to texture fractions (silt and clay) and OM (Table 3) by multiple regression. The PTFs were also tested differentiating between fine and coarse silt. The fine silt was significant (p < 0.05) in opposition to coarse silt (not significant). Although slightly better fitting models were obtained with the fine silt fraction, the regression on total silt was preferred since data on the silt fraction are more common than fine and coarse silt data.
Fig. 5a shows the relationship between the measured and PTF-predicted FC and WP. The General Model had good coefficients of correlation (FC = 72.3% and WP = 74.6%) but high RMSE (2.44 and 2.11). Fig. 5b, c, and d show the three submodels: Zone I had a better fit (in terms of R^; Table 3), especially for the WP, but Zones II and III showed weaker fits (in special Zone III due to the heterogeneous origin of its soils -valley bottoms and gypsiferous heights) than the General Model (lower R2adj). On the other hand, the RMSE for the submodels was reduced (in Zones I and II, not in Zone III). Thus, although the R2adj for the Distributed Model was lower than for the General Model, the
latter allows for FC and WP estimations with lower RMSE (in Zones I and II, where the main part of the irrigated crops are located) and is preferred to the General Model.
The concordance correlation coefficient was 0.85 and 0.89 for FC and WP respectively in the Global Model and 0.87 and 0.89 for the Distributed Model. Generally, concordance correlation coefficient lower than 0.90 are considered poor when applied to identifying biases in instrumental analysis (McBride, 2005); but for the kind of field data tested in this work, values close to 0.90 are deemed good enough.
4.2.2. Hydraulic Conductivity
The Saturated Hydraulic Conductivity in VID ranged from low to moderately low (Soil Survey Division Staff, 1993). Table 2 synthetize the results from the conductivity tests. The Ks data obtained showed a log-normal distribution. The t-test between the logarithmic transformations of Ks and Ks-Rosetta values (not provided) showed significant differences. Givi et al. (2004) explained that the soils (form subtropical climates) used to develop the software were quite different from typical arid or semi-arid soils, and Nguyen et al. (2015) pointed to the convenience of applying PTFs only for the regions where they were derived, which might explain the differences found. This fact evidences the need for accurate local PTFs.
Two kinds of PTFs were set up (Table 4): on Quantitative GC (through the normally distributed lnGC) and on Qualitative GC (dummy variable). The latter discriminated between high and low values of GC (4% was chosen as threshold; Fig. 2) and could be useful when quantitative GC data is not available but there are gypsum field tests available (e.g. precipitation field test with BaCl2). The PTFs were also linked to texture: sand or silt.
Table 2
Average results of the conductivity tests (Ks; m/d) performed by subgroup and particle size family. Different letters under the same column indicate significant differences (P < 0.05) between columns.
Particle size family
Fine textured soils
Coarse textured soils
Subgroup Horizon Fine (1) Fine Silty (2) Fine Loamy (3) Coarse Loamy (4) Loamy (shallow) (5) Loamy-skeletal (6)
Typic calcixerept A 0.118 0.159 0.617
B 0.068 0.183 0.513
Petrocalcic calcixerept A 0.140
B 0.113
Gypsic haploxerept A 0.107 0.516
B 0.193 0.167
C 0.176
Typic xerothent A 0.155 0.147 0.378 0.144
B 0.122 0.329 0.168
C 0.209
Typic xerofluvent A 0.079
B 0.078
ANOVA abc b cd d bcd a
Fig. 4. Homogeneous Zones for hydrological management in the Violada Irrigation District.
Table 3
Field capacity and wilting point (% weight) pedotransfer functions for the General Model, each homogeneous zone submodel and the Distributed Model; upon the variables clay, silt and organic matter (%); adjusted coefficient of determination (R^j), root mean square error (RMSE), and number of data (N) of each model.
Pedotransfer functions R2dj (%) RMSE N
General Model (whole VID)
FCvid — 5.701 + 0.211 * Clay + 0.234 * Silt 72.3 2.44 88
+ 1.344 * OM
WPviD — - 2.536 + 0.269 * Clay + 0.149 * 74.6 2.11 88
Silt + 1.212 * OM
Distributed Model (I + II + III)
Submodels Zone I FCi = 7.506 + 0.264 * Clay + 75.1 1.96 27
0.139 * Silt + 1.341 * OM
WPI = - 2.705 + 0.327 * Clay + 86.5 1.45 27
0.105 * Silt + 1.746 * OM
Zone II FCII = 21.665 + 0.1297 * Clay + 54.7 1.71 25
1.573 * OM
WPII = 8.807 + 0.163 * Clay + 57.8 1.93 25
1.761 * OM
Zone III FCIII = 11.608 + 0.207 * Clay + 33.2 2.71 36
0.173 * Silt
WPIII = - 2.099 + 0.294 * Clay 53.4 2.37 36
+ 0.158 * Silt
Joint Distributed Model
FCi + ii + Hi 54.2 - 88
WP, + ii + iii 66.9 - 88
The Ks PTFs for the VID were disaggregated by the three homogeneous zones for fine textures (Table 4). Loamy (shallow) and Loamy-skeletal textures (Fig. 4) were fitted with a single Ks model since these textures were not linked to GC and the coarse fragments (not considered for the PTF) affect the soil hydraulic properties (Rawls and Brakensiek, 1989, Poesen and Lavee, 1994).
Fig. 6 shows the relationship between the measured and PTFs-predicted LnKs values, together with the residual distribution and residual histograms for each model. The General Model (whole VID) presented R2adj of 36.9% for the Quantitative GC and 34.1% for the Qualitative GC and a RMSE of 0.76 and 0.75 respectively (Table 4). The residuals were normally distributed but with a wide histogram. The Distributed Model increased the R2adj to 53.8% and 47.9% for the Quantitative and Qualitative GC PTFs respectively. The worst adjustments were obtained for Zone II with R2adj of 23.8% and 20.0% for the Quantitative and Qualitative GC models respectively. The low variability of the only significant variable for Zone II (lnGC) may explain the low R2 coefficient. On the other hand, for the other fine textured soils the models presented better R2adj (Quantitative GC: 68.0% in Zone I and 69.1% in Zone II; Qualitative GC: 53.4% and 68.0%) than the General Model. The RMSE of the distributed model were also improved for Zone I (0.58 in the Quantitative GC model and 0.64 in the Qualitative GC model) and Zone III (0.51 and 0.52); but not in Zone II (0.90 and 0.93). The RMSE was especially low for the Loamy (shallow) and Loamy-skeletal textures (0.15) (Table 4).
Fig. 5. Predicted and measured values of FC and WP for General Model (a) and the submodels for the distributed model (b, c and d).
The use of the Distributed Model instead of the General Model allows for improving the fits, as shown by the R^j (36.9% for Quantitative GC and 34.1% for Qualitative GC to 53.8% and 47.9% respectively) and especially the distribution of the residuals (Fig. 6). The Distributed Model is thus preferred to the General Model in the estimation of Ks.
4.3. Salinity and Sodicity
Sulfates were the main source of salts in VID (saturated paste extracts) with a strong relation (slope ~ 1) between the estimated SO42 - and the sum of Ca2+ and Mg2+ pointing to a common origin in sulfate salts. The EC-Cluster showed three clusters (from the analysis of the principal components on the variables ECe and EC1:5; Fig. 7a). The
first cluster (named Non-gypsum Saturated) represents the horizons whose saturated paste extract was not saturated in gypsum. The second cluster (Gypsum Saturated) represents the horizons saturated in gypsum. In both clusters, the SARe remained very low, as shown by the isolines calculated from the multiple regression of SARe on the ECe and EC1:5 (Fig. 7a).The third cluster (Saline) corresponds to the horizons of the three saline pits (C-29, C-32 and C-33; Fig.2) and presented higher SARe. In the 71 samples (N) from the first and third cluster, the ECe was linearly related to EC1:5:
ECe = 0.37 + 4.68 X ECVi — ; R2 = 86.3%; RMSE =1.21—; N =71
These coefficients were similar to those obtained by Herrero and Bercero (1991). But for the samples in the Gypsum Saturated cluster, the ECe showed no relationship to EC1:5 (coefficient of regression not
Table 4
Pedotransfer functions for the natural logarithm of the hydraulic conductivity in m/d (LnKs) on the variables sand, silt and GC (%): adjusted coefficient of determination (R2adj), root mean square error (RMSE), and number of points used (N) for the general model, submodels (by zones and textures) and distributed model.
Quantitative gypsum content
Qualitative gypsum content
N Pedotransfer function R^di (%) RMSE GC (%) Pedotransfer function R^di (%) RMSE
General Model (whole VID) 69 LnKs = - 2.850 + 0.033 * Sand + 0.385 * LnGC 36.9 0.76 < 4 LnKs = - 2.695 + 0.026 * Sand 34.1 0.75
> 4 LnKs = - 1.702 + 0.026 *
Distributed Model (I + II + III)
Submodels Fine Zone I 17 LnKs = - 1.704 - 1.576 * LnGC 68.0 0.58 LnKs = - 2.759 + 0.035 * Sand 53.4 0.64
Zone II 16 LnKs = - 2.547 + 0.442 * LnGC 23.8 0.90 <4 LnKs = - 2.471 20.0 0.93
>4 LnKs = - 1.429
Zone III 23 LnKs = - 0.018 * Sand - 0.033 * Silt + 0.824 * 69.1 0.51 <4 LnKs = - 1.136 - 0.022 * Silt 68.0 0.52
LnGC >4 LnKs = 1.055 - 0.022 * Silt
Coarse 13 LnKs = - 1.535 - 0.011 * Sand 48.9 0.15 LnKs = - 1.535 - 0.011 * Sand 48.9 0.15
Joint Distributed Model 69 53.8 - 47.9 -
Fig. 6. Predicted and measured values of Ln(Ks) for the general and distributed models with Qualitative and Quantitative GC. The graphs on the right show the residuals distribution and histogram for each model.
Fig. 7. Relationships between electrical conductivities; SARe and GC for the clusters obtained from the principal components analysis for all pit horizon samples: a) relationship between ECi:5 and ECe with the isolines calculated for the SARe; b) relationship between the GC and the ECi:5.
significant):
EC, = 2.96 — ;
RMSE = 0.43—; N = 39
Fig. 7b shows the effect of the GC on EC1:5. The electrical conductivity of the solution 1:5 was not affected by the GC in the Non-Gypsum Saturated cluster: the EC1:5 was below 1 dS/m and the GC below 4%. In the Gypsum Saturated cluster the EC1:5 increased with the GC for values below 4% and over this threshold, EC1:5 remained constant (~2.2 dS/m) regardless of the GC. Similar results were obtained in close-by irrigated areas by other authors (Abrisqueta et al., 1962; Herrero et al., 2009; Moret-Fernandez and Herrero, 2015; Nogues et al., 2006).
The Saline cluster was found on two locations: (i) on the gypsum heights in the south (C-32 and C-33; Fig. 1) and (ii) to the Southeast of the glacis (C-29; Fig. 1). Pits C-32 and C-33 showed high GC (~12%), but the extract was not saturated in gypsum due to the presence of more soluble salts as NaCl (182 Cl- mmolc/l and 157 Na+ mmolc/l at the A horizon in C-32 pit) rising the gypsum solubility (Casby-Horton et al., 2015). Only C-32 pit points to a potential problem of salinity (ECe = 23dS/m) or sodicity [SARe = 22.6 (mmolc/l)0 5] but the high
Table 5
Mean, standard deviation (SD) and number of samples (N) of ECe (dS/m) and SARe [(mmolc/l)0 5] by subgroups and horizons.
ECe (dS/m) SARe (mmolc/l)0-5
Subgroup Horizon Mean SD Mean SD N
Typic & Petrocalcic A 1.90 1.22 1.02 0.89 14
Calcixerept
B 2.31 1.65 1.65 1.44 13
C 3.03 1.78 1.59 2.05 4
Gypsic Haploxerept A 2.27 0.93 0.60 0.26 9
B 2.58 0.73 0.65 0.27 17
C 2.82 0.26 0.49 0.24 7
Typic Xerorthent A 4.38 6.18 2.76 6.38 12
B 4.87 4.52 4.63 7.75 9
C 2.77 2.60 - 2.57 6
Typic Xerofluvent A 1.70 0.86 0.63 0.29 5
B 2.33 0.95 0.81 0.29 10
C 3.58 - 1.10 - 1
Table 6
Correlation coefficients between the salinity variables and factors. Coefficients higher than 0.75 are typed bold-face.
Salinity factor Carbonates-Gypsum factor Bicarbonate factor
CCE - 0.13 - 0.79 0.23
GC - 0.05 0.86 0.15
ECe 0.97 0.16 0.13
Cl- 0.93 - 0.03 0.00
HCO3- - 0.09 0.03 - 0.92
Na+ 0.96 - 0.01 0.05
Ca2 + 0.35 0.55 0.46
Mg2 + 0.83 0.25 0.18
Ca2+ was not related to any factor in particular, with similar correlation coefficients with the three factors, likely due to the high presence of this anion all over VID (Fig. 3).
Fig. 8 presents the horizon samples on the plane of the first two factors classified by the EC-cluster analysis. The three clusters are much differentiated in the graph as expected. The Saline Cluster spread widely along the Saline Factor, being markedly different from the other samples in this factor. The other two clusters spread along the Carbonates-Gypsum factor axis with positive values for the Gypsum Saturated cluster samples and negative values for the Non-Gypsum Saturated cluster samples.
4.4. Linking hydrologic features to irrigation practice and model applications
position and the absence of water logging problems close to C-32 or C-33 points to the original material as the primary cause of the salinity in this area. On the contrary, the slightly saline pit C-29 (Southeast of the glacis) had low GC (1.8%) and ECe > 4 dS/m in the lower horizons. In this case, the location in a relative low-lying position points to poor drainage as responsible of the salinity.
On the other hand, the areas prone to water logging (Gypsic Haploxerept with imperfect drainage to the Northwest of VID, Fig. 1) did not show evidences of salinity (ECe < 4 dS/m). This together with the average ECe and SARe by subgroup and horizon (Table 5), made it clear that there is currently no generalized problem of salinity at the VID associated to lack of drainage, probably due to the drainage network implemented by the farmers since irrigation was established. Only the Typic Xerorthent subgroup shows the mean ECe above the salinity limit (4 dS/m) in the A and B horizons, due to the ECe in pit C-32 [if C-32 is excluded the average ECe's for the subgroup were: 2.9 dS/ m (A horizon); 2.58 dS/m (B); and 2.77 dS/m (C)].
The variability of the salinity in VID was explained by Factor Analysis using as variables CCE, GC, ECe, and Na+, Ca2+, Mg2+, Cl-, and HCO3- concentrations in the saturation extract. The analysis explained 81.3% of the variance with three independent factors linked to (1) ECe and cations, (2) CCE and GC and (3) HCO3-. The first factor was labelled as Salinity Factor and accounts for 48.2% of the variance, positively correlated to ECe,Cl + ,Na + and Mg2+ (Table 6), showing the strong relationship between the ECe and the presence of soluble salts as NaCl; and the independence of salinity from the GC. The second factor accounted for 20.1% of the variance and was linked negatively to CCE and positively to GC; and labelled Carbonates-Gypsum Factor (a factor also found in the companion paper). The third factor was only linked negatively to HCO3-; therefore, it was labelled as Bicarbonate Factor.
Fig. 8. Factor scores of each horizon identified by their cluster membership (Non-Gypsum Saturated, Gypsum saturated and Saline) in the two first factors. The circle at the top-right corner shows the coefficients of correlation of the variables used in the analysis with the two first factors.
The soil hydrological features at VID define three hydrological homogenous zones in the VID: Zone I corresponds to glacis (with the highest CaCO3 accumulation), Zone II to valley bottoms high in gypsum, and Zone III to valley bottoms without gypsum, as well as colluvial slopes and some hills high in GC. Zone I was characterized by lower AWC and higher Ks in relation to Zones II and III. Plots in Zone I are expected to have higher infiltration capacity than Zone II and III, a point to be considered during the irrigation management decision process.
Those zones allow setting up two models in the VID and a collection of PTFs for the soil hydrological features (FC, WP and Ks). These two models are: (i) General Model (aggregated; only one PTF per feature for the whole VID) and (ii) Distributed Model (sum of the submodels for each zone) with PTFs defined for the three zones. The Distributed Model may be used in modeling nearby areas of similar characteristics to predict FC, WP or Ks from texture and GC providing reliable estimates. Both the RMSE and the distribution of the residuals (in the Ks) suggest that using Distributed Models may yield better results than the General Model (lower RMSE for the main part of the surface) in establishing hydrologic properties. And these estimated hydrologic properties may help devising irrigation management strategies (doses and frequencies) suited to the actual soil properties. The differences in Ks estimates from a generalist PTF model (like Rosetta) and the PTFs obtained for the VID evidence the need to develop local PTFs for Ks estimation and to take into account other, less general variables than used in generalist models. These locally established PTFs should be based on the actual soil type distribution in the area, especially if the result will be the basis to develop water movement models. And clearly gypsum should be considered in developing PTFs for Ks in arid or semiarid areas.
The areas with lower Ks (also with heavier textures and low lying position in the landscape) could be more prone to salinity development induced by irrigation (although only one pit in these areas has shown some salinity, C-29). Thus, the drainage network of the lowlands should be maintained and the attention should be paid to potential salinity development in these areas. The high levels of Ca2+ in the soil solution (derived from the prevalence of gypsum in the area) yield quite low SARe in most soils (slightly higher in the Typic Xerorthents, due only to saline pits C-32 and C-33). Soil infiltration problems may arise anyway from the use of high quality waters (low salinity) and sprinkler irrigation that could lead to the formation of soil crust. However, this problem could be skipped by keeping higher water content in the upper soil through more frequent irrigations.
5. Conclusions
This work completes the analysis started with the making of the VID Soil Map and shows the importance of the GC in any further analysis to be developed in the VID (use of soil water models e.g.) or nearby areas with similar characteristics.
The differences in hydrologic properties down to 60 cm (root depth
observed in the field) were used to define three homogeneous zones in regard to water storage capacity; and one additional zone (coarse textured soils) with special characteristics in terms of Ks. The calcic glacis (Zone I) with coarser textures may need more frequent irrigation, due to their lower AWC than the valleys bottoms (Zone II and III only differentiated by their GC) with heavier textures, lower position in the landscape, and limited permeability which are more prone to saliniza-tion.
The zones were used to establish specific PTFs for FC and WP (% weight) and Ks in each one The PTFs reveal that FC and WP (in % weight) depend on texture and OM while point to the relevance of GC on the Ks. In shallow soils or skeletal textures, GC had no effect on Ks. The Distributed Models (PTFs for each zone) were preferred to the General Model (PTFs for the whole VID) as they resulted in lower RMSE of the residuals for most of the VID area. These results also point to the need to establish local PTFs, disaggregated for soil units of similar characteristics when possible; and suggest the use of gypsum (quantitative or quantitative) as a variable in the PTFs for Ks in arid, gypsum-rich environments.
The analysis of the ECe and SARe (considering the GC) suggested that there are no significant salinity or sodicitiy problems, at the moment, possibly due to the artificial drainage network implemented by the farmers since the 1940s. At present, salinity in VID seems to be linked to the original materials rather than to limited drainage. Nevertheless salinity is a latent risk in the area that should be always considered in irrigation planning to ensure adequate salt leaching, particularly in valley bottoms with low Ks.
Sulfates were the main source of dissolved ions in the soils of VID and the saturated extract is dominated by calcium. The relationship between ECe and EC1:5 was completely different in the horizons with (ECe constant) and without gypsum (linear increase with EC1:5). The high levels of Ca2 + (and Mg2+) in the soil solution and the presence of gypsum may prevent future sodicity issues.
In future works, these PTF models, once validated, will provide a tool for assessing the environmental effect (on return flows and soil properties) of irrigation and crop management or climate change scenarios in VID or in in nearby areas with the same type of soils and low detail soil data (e.g. potential modernization districts or new irrigated areas). Also, the low permeability and silty texture of many valley soils along with the sprinkler irrigation with high quality (low salinity) waters calls for the preventive monitoring of infiltration.
Acknowledgements
This work was funded by the Spanish Ministry of Science and Innovation through the research grants AGL2010-21681-C03-03 and AGL2013-48728-C2-2-R and by the Ebro River Basin Authority through a collaboration protocol. Thanks are given to the Almudevar Water User Association and the Ebro River Basin Administration for their support.
Appendix A. Supplementary data
Supplementary data associated with this article can be found in the online version, at http://dx.doi.org/10.10167j.geoderma.2017.04.024. These data include Google map of the most important areas described in this article.
References
Abrisqueta, C., Guillén, M.G., Fernández, F.G., Caro, M., 1962. Contribución al estudio de la determinación de la salinidad del suelo. Anales de Edafología y Agrobiología XXI 545-554.
Aragüés, R., Millán, M., Quilez, D., Fernández, M., 1986. Métodos de medida de la
salinidad del suelo. I y II. M.A.P.A. Insitituto Nacional de Investigaciones Agrarias, Madrid.
Artieda, O., Herrero, J., Drohan, P.J., 2006. Refinement of the differential water loss method for gypsum determination in soils. Soil Sci. Soc. Am. J. 70, 1932-1935.
Bouma, J., 1989. Using soil survey data for quantitative land evaluation. In: Stewart, B.A. (Ed.), Advances in Soil Science. Springer US, New York, NY, pp. 177-213.
Casby-Horton, S., Herrero, J., Rolong, N.A., 2015. Chapter four - gypsum soils—their morphology, classification, function, and landscapes. In: Donald, L.S. (Ed.), Advances in Agronomy. Academic Press, pp. 231-290.
CBDSA, - Comisión del Banco de Datos de Suelos y Aguas, 1983. SINEDARES, Manual para la Descripción Codificada de Suelos en el Campo. MAPA, Madrid, pp. 137.
Díaz, L., Herrero, J., 1992. Salinity estimates in irrigated soils using electromagnetic induction. Soil Sci. 154, 151-157.
Franzmeier, D.P., 1991. Estimation of hydraulic conductivity from effective porosity data for some Indiana soils. Soil Sci. Soc. Am. J. 55, 1801-1803.
Givi, J., Prasher, S.O., Patel, R.M., 2004. Evaluation of pedotransfer functions in
predicting the soil water contents at field capacity and wilting point. Agric. Water Manag. 70, 83-96.
Harman, H.H., 1967. Modern Factor Analysis, second ed. Univ. of Chicago Press, Oxford, England, pp. 474.
Herrero, J., Artieda, O., Hudnall, W.H., 2009. Gypsum, a tricky material. Soil Sci. Soc. Am. J. 73, 1757-1763.
Herrero, J., Bercero, A., 1991. La salinidad en el nuevo regadío de Quinto (Zaragoza). Suelo y planta 602, 585-602.
Ibarra, P., 2004. La diversidad edáfica del territorio aragonés. In: Peña Monné, J.L., Longares Aladrén, L.A., Sánchez Fabre, M. (Eds.), Geografía física de Aragón. Aspectos generales y temáticos. Institución Fernando el Católico and Universidad de Zaragoza, Zaragoza, pp. 41-53.
Isidoro, D., 1999. Impacto del regadío sobre la calidad de las aguas superficiales del Barranco de La Violada: salinidad y nitratos. PhD Thesis, Univerisidad de Lleida, Lleida. (267 pp).
Jiménez-Aguirre, M.T., Isidoro, D., Usón, A., 2017. Soil variability in La Violada Irrigation District (Spain): I Delineating soil units for irrigation. Geoderma. http://dx.doi.org/ 10.1016/j.geoderma.2017.04.025. (in press).
Lin, L.I.-K., 1989. A concordance correlation coefficient to evaluate reproducibility. Biometrics 255-268.
MARM, 2002. Ministerio de Medio Ambiente Medio Rural y Marino. In: Plan Nacional de Regadíos.
MARM, 2006. Ministerio de Medio Ambiente Medio Rural y Marino. In: Plan de choque de modernización de regadíos.
McBride, G., 2005. A proposal for strength-of-agreement criteria for Lin's concordance correlation coefficient. In: NIWA Client Report: HAM2005-062.
Minasny, B., Hartemink, A.E., 2011. Predicting soil properties in the tropics. Earth Sci. Rev. 106, 52-62.
Moret-Fernández, D., Herrero, J., 2015. Effect of gypsum content on soil water retention. J. Hydrol. 528, 122-126.
Nguyen, P.M., Van Le, K., Botula, Y.-D., Cornelis, W.M., 2015. Evaluation of soil water retention pedotransfer functions for Vietnamese Mekong Delta soils. Agric. Water Manag. 158, 126-138.
Nogués, J., Herrero, J., 2003. The impact of transition from flood to sprinkler irrigation on water district consumption. J. Hydrol. 276, 37-52.
Nogués, J., Robinson, D.A., Herrero, J., 2006. Incorporating electromagnetic induction methods into regional soil salinity survey of irrigation districts. Soil Sci. Soc. Am. J. 70, 2075-2085.
Oosterbaan, R., Nijland, H., 1994. Determining the saturated hydraulic conductivity. In: Ritzema, H. (Ed.), Drainage Principles and Applications. Ed. 2 ILRI Publication 16, Wageningen, The Netherlands, pp. 435-476.
Pachepsky, Y.A., Rawls, W.J., Lin, H.S., 2006. Hydropedology and pedotransfer functions. Geoderma 131, 308-316.
Poesen, J., Lavee, H., 1994. Rock fragments in top soils: significance and processes. Catena 23, 1-28.
Rawls, W.J., Brakensiek, D.L., 1989. Estimation of soil water retention and hydraulic properties. In: Morel-Seytoux, H.J. (Ed.), Unsaturated Flow in Hydrologic Modeling. Springer, Netherlands, pp. 275-300.
Schaap, M.G., Leij, F.J., van Genuchten, M.T., 2001. ROSETTA: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol. 251, 163-176.
Simunek, J., Sejna, M., Saito, H., Sakai, M., van Genuchten, M.T., 2008. The Hydrus-1D software package for simulating the movement of water, heat, and multiple solutes in variably saturated media, version 4.0, HYDRUS software series 3. In: Department of Environmental Sciences, University of California Riverside, Riverside, California, USA.
Slatni, A., 1996. Elaboration et evaluation des alternatives pour l'amelioration de
l'utilisation de l'eau au sein de la communaute d'irrigants d'Almudévar. Master Thesis, Instituto Agronómico Mediterráneo de Zaragoza (CIHEAM-IAMZ), Zaragoza (132 pp).
Soil Survey Division Staff, 1993. Soil survey manual. In: Soil Conservation Service. United States Department of Agriculture Handbook, pp. 18.
Soil Survey Staff, 2011. Soil Survey Laboratory Information Manual. Soil Survey
Investigations Report No. 45, Version 2.0.R. Burt (ed.). United States Department of Agriculture. Natural Resources Conservation Service.
Soil Survey Staff, 2014. Keys to Soil Taxonomy, 12th ed. United States Department of
Agriculture. Natural Resources Conservation Service. Stambouli, T., 2012. Gestión avanzada del riego por aspersión en Parcela: Aplicación en el
Valle Medio del Ebro. Universidad de Zaragoza, Zaragoza, Master Thesis (190 pp). Szabolcs, I., Várallyay, G., 1979. Review of Research on Salt-Affected Soils. Unescopp. 137.
USSLS, United States Salinity Laboratory Staff, 1954. Diagnosis and improvement of
saline and alkali soils. In: Agriculture Handbook No 60, pp. 160. Wagner, B., Tarnawski, V.R., Hennings, V., Müller, U., Wessolek, G., Plagge, R., 2001. Evaluation of pedo-transfer functions for unsaturated soil hydraulic conductivity
using an independent data set. Geoderma 102, 275-297. Wessolek, G., Bohne, K., Duijnisveld, W., Trinks, S., 2011. Development of hydro-pedotransfer functions to predict capillary rise and actual evapotranspiration for grassland sites. J. Hydrol. 400, 429-437. Wösten, J.H.M., Lilly, A., Nemes, A., Le Bas, C., 1999. Development and use of a database
of hydraulic properties of European soils. Geoderma 90, 169-185. Wösten, J.H.M., Pachepsky, Y.A., Rawls, W.J., 2001. Pedotransfer functions: bridging the gap between available basic soil data and missing soil hydraulic characteristics. J. Hydrol. 251, 123-150.
