Agricultural practices and geostatistical evaluation of nitrate pollution of groundwater in the Júcar River Basin District , Spain

doi: 10.9755/ejfa.2015-04-346


INTRODUCTION
The Júcar River Basin District is located in the Communities of Castilla La Mancha and Valencia, Spain (CHJ 2007, 2009, 2011;Ferrer et al., 2012) (Fig. 1).It extends over about 42 851 km 2 in a hilly area with elevations ranging from sea level to 2000 m.Intensive agriculture is an important land use activity, occupying about 20735 km 2 of the cultivated area, of which 15% is irrigated.Due to adverse climatic conditions, water resources are of considerable importance.The largest portion of the basin has a semi-arid climate, which reduces the availability of surface water.This necessitates careful overall management of water sources to ensure effective use of water by all users, including urban consumers.Tourism and industry are large water consumers, but agriculture is the most significant.On the basis of data from 2005, consumption of water for agricultural purposes is estimated at 2452 hm 3 in the Júcar River Basin District.Beside this large consumption, there are significant problems with the quality of groundwater, including nitrate pollution, especially close to the most intensive agriculture areas.
For this study, a space-time table of nitrate concentrations obtained from 6381 monitoring locations within the Júcar River Basin District between 1969 and 2012 was used.Nitrate concentrations reach 300 mg/l in several locations, compared to the 50 mg/l maximum recommended concentration in water for human consumption (Nitrates Directive 91/676/EEC; Council of the European Communities, 1991).The major source of nitrates is fertilizers used in intensive agriculture, and as nitrates are highly soluble, they are leached into water-courses and lakes, creating serious pollution problems in surface water and groundwater (Gill, 1996;Guerrero, 2003;Mejías et al., 2012).
The type of soil determines the extent to which nitrates are leached, by influencing the surface flow and infiltration rate (Brito et al., 2006;Lencastre and Franco, 2006;Silva and Almeida, 2009).More permeable formations like sand, gravel and karst are more likely to conduct large amounts of water and leach pollutants such as nitrates.In general, water flows through the different layers of topsoils when rainfall exceeds evapotranspiration and soil moisture, which normally occurs when vegetation cover is low or nonexistent.

GEOGRAPHICAL AND GEOLOGICAL FRAMEWORK OF THE JÚCAR RIVER BASIN DISTRICT
The study area comprises the hydrographic basins and sub-basins draining into the Mediterranean Sea from the left bank of the mouth of the river Segura to the mouth of the river Sénia including the basin of the latter, and the endorheic basin of Pozohondo (Fig. 1).Of the total area of the District (42851 km 2 ), 50% is occupied by seminatural forest and 46% by agriculture, the latter playing an important role in the environmental conditions of the region.
Both irrigated and non-irrigated crop types are found in the Júcar River Basin District.In irrigated crop types, the production of vegetables and fruit predominates, and water is supplied by surface sources and aquifers.In non-irrigated crop types, in which cereals predominate, rainwater is needed.The fertilizers used for both crop types are identical and nitrate-based.Due to the difference in water supply and needs, irrigated agriculture plays a major part in nitrate leaching to the soil, and thus makes a more significant contribution to groundwater contamination.In non-irrigated crop types, despite the use of similar nitrate-based fertilizers, nitrates are only leached by rainfall (Suarez-Varela et al., 1998).
From a geomorphological point of view, the region has three major structures: mountain systems, the continental plateau and the coastal plain (Fig. 2).Of the mountain systems, the most important are the Iberian System in the north, reaching 2000 m, and the Betic System in the south and southeast, with elevations around 1200 m.The continental plateau west of the Betic and Iberian systems is a plain with elevations around 700 m.Finally, the coastal plain has mainly alluvial soils and is delimited to the north by the Iberian System, to the south by the Betic System, to the west by the continental plateau and to the east by the Mediterranean Sea.
Geologically, carbonate rocks prevail, mostly calcarenites and marls, especially in the mountains.In a few sectors other lithologies are also found, mainly detrital, siliceous detrital, selenitic and clayish.All the rocks are sedimentary of the Mesozoic age.Occasionally, some Paleozoic outcrops and some materials of the Tertiary age occur, such as alluvial deposits on river banks.Groundwater is of primary importance given that 73% of all water used in the region, i.e., about 2268 hm 3 /yr, is provided by aquifers (Custodio and Llamas, 1996).Although they are supplied from various sources, some aquifers have already been overexploited (CHJ, 2007).

METHODS
The data table used in this case study links several networks of wells and boreholes/piezometers built or maintained for water exploitation and monitoring purposes.There were 6381 monitoring stations available for study located in 113 groundwater bodies, totalling 24276 records (Fig. 3).The time series records begin in 1969 and continue to 2012 (41 hydrological years), and during this period each monitoring location recorded nitrate concentrations together with other parameters.Due to alternating activity and inactivity of the monitoring stations, the time series records are not complete for any of the monitoring locations.
Before the data were analysed, groundwater bodies with low permeability or of only local interest were excluded from the study (Fig. 4).Then, for the remaining groundwater bodies with higher permeability and therefore exploitation capacity, the five most important monitoring locations were selected, defined as those with most temporal and recent records.Using these criteria, 275 monitoring locations were selected out of a total of 6381 available.Fig. 5 displays the number of measurements of nitrates for each hydrological year used in the case study.As can be observed, the number of yearly observations is highly variable, ranging between a few records to almost 400 in the hydrological year 1999/2000.
The first statistical analysis consisted of evaluating the linear correlation coefficient () to the representation of nitrate concentrations vs. sampling date.The signal of  (positive, negative or null) represents the level and type of correlation.A  close to zero means that no temporal trend exists, a positive one means that an increase in nitrates over time is observed and a negative  means an inverse trend.A higher absolute value of , close to unity, means a consistent increase or decrease in concentrations; a lower  means a more heterogeneous and less consistent variation over time.
Based on the results of this analysis, balancing the quantity of sampling points and variations in nitrate concentrations over time, each groundwater body was classified as one of the following: Sharp Rise, Gentle Rise, Inconclusive, Gentle Fall and Sharp Fall.The classes Sharp Rise and Sharp Fall were adopted for  higher than 0.7, meaning greater robustness of the trends.The intermediate classes Gentle Rise and Gentle Fall include data for which  ranges between 0.7 and 0.5.The remaining class of Inconclusive was applied to  results lower than 0.5.The sign of the slope distinguishes between Rise (positive sign) and Fall (negative sign) classes.This analysis shows that variations in nitrate concentrations over time for each groundwater body are usually consistent between different sampling points, which can be explained by more regional sources of variations in the dependent parameters.
Possible regional sources may be different rainfall patterns, diversity in land uses, diversity in fertilizers used in agriculture, expansion of crops and variations in the quantity of water extracted from wells and boreholes.The data available for this study only enabled analysis of the relationship between temporal variations in nitrate concentrations and yearly rainfall for each groundwater body.Thus, a bivariate statistical analysis between total yearly rainfall and nitrate concentrations for each groundwater body was performed and the results were characterized by a linear correlation coefficient.Leaching and transport of nitrates from the surface to groundwater bodies takes a certain time, which depends on rainfall and the types of soil and rock involved.
To take into account delays between rainfall and transport of nitrates through the rocks of the mainly detrital and carbonate geological formations with different hydraulic conductivities, correlations were analysed between the observed concentrations and rainfall in previous years up to a maximum of two years.The magnitude and the sign of the linear correlation coefficient determine the type and level of dependency of these two variables.If the coefficient is positive and close to one, it means a high and positive correlation, so the nitrate concentration increases as rainfall also increases.On the other hand, if the sign is negative, there is an opposite correlation between nitrate concentration and rainfall.Low coefficient values (close to zero) mean that the variables of nitrate concentration and rainfall are independent.
The infiltration rate of water and the time taken to reach groundwater bodies depends on the hydraulic conductivity of the rocks.In detrital rocks, water percolation follows the small and interconnected pores of the rock and therefore capillarity is an important factor affecting flow.By contrast, in carbonate rocks the percolation rate is higher than in detrital rocks because of its secondary porosity and permeability caused by fractures and dissolution of carbonate components that create sinkholes and caves.The geological nature of the different materials that compose each groundwater body is crucial for interpretation of the results.Therefore, the set of groundwater bodies were grouped according to their geological characteristics, as To optimise the analysis of the correlation between nitrate concentrations and climate (rainfall), the annual recharge to each groundwater body in the Júcar River Basin as calculated by the PATRICAL simulation model was also used (Pérez-Martín, 2005).The PATRICAL model quantifies the amount of recharge by taking into account a set of contextual parameters in the calculations, giving more realistic results for the amounts of water that reach the ground and infiltrate into the lower layers.For example, in a rainy year, the usual yearly estimation does not take monthly variations in rainfall into account, but PATRICAL calculates the monthly recharge, the annual recharge being the sum of the monthly amounts.The recharge is computed as rainfall minus surface flow and evapotranspiration.The PATRICAL model enables more effective management and planning of the water resources in the region.
After this statistical analysis, which assessed possible relationships between nitrate concentrations and some regional contextual variables (precipitation and type of groundwater bodies), the study proceeded with the construction of annual average images of nitrate concentrations in the area of study (Júcar River Basin).To create average yearly images, a geostatistical approach was used, which included (Soares, 2000;Paralta and Ribeiro, 2003;Guimarães et al., 2014;Quental et al., 2012;Roxo et al., 2016): (a) calculation and fitting of variograms of nitrate concentrations in space (2D) and time (third dimension); (b) generation of various simulated scenarios of nitrate concentrations for each hydrological year using Sequential Gaussian Simulation (SGS); (c) determination of the average of the simulated nitrate concentrations; (d) evaluation of the areas above and below the threshold of 50 mg/l nitrates.The space-time approach is effectively a 3D model, meaning that the X and Y coordinates represent 2D space and time represents the third dimension (usually Z).For the estimation of each location the nearest data, in space and in time, were used (Almeida et al., 2004).
For geostatistical modelling purposes the geoMS software and a parallelised version of the sequential Gaussian simulation were used (Alexandre and Almeida, 1998;Nunes and Almeida, 2010).Due to the broad dispersion of data with no preferential orientation, it was decided to compute averages for the variograms encompassing all directions together (omnidirectional variograms) (Isaaks and Srivastava, 1989;Goovaerts, 1997).Also, to estimate yearly average images, simulations plus average of simulations were performed instead of kriging for declustering purposes, thus avoiding the artefacts in the images that would have been obtained by kriging (Almeida and Lopes, 2005).For processing simulated images it is also better to provide a binary image above and below a threshold, which is done at the end of the study (Almeida, 2010).

Correlation analysis between nitrate concentrations and the contextual variables rainfall, infiltration and groundwater body major rock type
The temporal trends in nitrate concentrations by groundwater body are displayed in Fig. 7.The 24 276 records of nitrate concentrations from samples collected in 6 381 locations (wells, boreholes and piezometers) located in 113 groundwater bodies show irregular trends of temporal variations, some clearly positive or negative, but most inconclusive.This trend analysis was based on the best fit slope of trend lines, and groundwater bodies were classified as Rise or Fall, Sharp (for  above 0.7), Gentle (with  between 0.7 and 0.5) or Inconclusive (for values of  below 0.5).The agreement between concentrations over time in various locations could be due to any regional cause affecting the groundwater bodies similarly, for instance the occurrence of dry or wet periods, or undetermined causes.
Trends in concentrations are undoubtedly related to land use, as can be observed in the continental plateau and the coastal plain, in both of which areas intensive agriculture is the major land use.
Analysis of the correlation between rainfall and nitrate concentrations began with bivariate diagrams between these variables for each groundwater body and for each hydrological year.Average nitrate concentrations for each groundwater body and for each hydrological year were computed for this analysis.
The highest correlation coefficient was observed for carbonate crystalline rocks, but did not exceed -0.50.If the correlation coefficients are calculated with the recharge values provided by PATRICAL model, slightly better correlations are observed.Even more, the calculations of correlations between nitrate levels in a specific year and the yearly rainfall of previous years (one and two previous years) did not yield substantially higher correlation coefficients.Even so, for carbonate crystalline groundwater bodies, negative correlations prevailed, implying a dependence between nitrate dispersion rate and lithology.In carbonates, water flows preferentially along fractures in which the flow velocity is much higher, resulting in more effective dispersion of pollutants.
Regarding carbonate detrital-type rocks, the correlation between nitrate concentrations and recharge is slightly better than that obtained for rainfall, but is far from good.The use of time delays in this bivariate analysis improves correlation coefficients in some cases.As was found for carbonate crystalline rocks, more groundwater bodies are observed with negative correlations than with positive correlations.Correlations with recharge values instead of rainfall show small, but not significantly greater improvements in results.The reasons can be assumed to be the same as those identified above for carbonate crystalline bodies: in carbonate rocks flow rates are high, resulting in dispersion of pollutants over a larger volume of water and more quickly.
Finally, in detrital rocks there are more groundwater bodies with positive correlation coefficients than with negative ones.A positive correlation coefficient means that an increase in seepage results in an increase in nitrate concentrations reaching the groundwater body.This may be a consequence of the hydraulic properties of detrital rocks, in which water flows between connected pores with lower conductivity rates than in carbonate rocks.
If infiltration and flow rates are lower, the increases in nitrate concentrations are more local, and for this reason conflicting results are obtained between close sampling points for the same groundwater body.
Table 1 displays the positive and negative correlations observed by groundwater body lithology type.To summarise, the linear correlations obtained between yearly rainfall or recharge calculated by the PATRICAL model and nitrate concentrations show mainly inverse

The geostatistical space-time model of nitrate concentrations
The geostatistical approach begins with the calculation of space and time variograms of nitrate concentrations and is followed by the stochastic simulation of nitrate concentrations using SGS.
As stated above, spatial variograms were computed merging all directions together (omnidirectional variograms); time variograms represent the variations in nitrate concentrations over time for each location.The spatial variograms show an effective correlation distance or range up to a maximum of 50 km; for the temporal calculations, a 500 year range of correlation was computed, mostly based on extrapolation of the results of the variograms.Fig. 8 illustrates the experimental variograms and the theoretical models fitted.
For the SGS, 350 realizations were made for each hydrological year, and the average of each set of 350 images was computed.Typically, each realization of a variable is a possible scenario of the studied variable, and honours the data in the measurement locations and the univariate statistics and variograms of the true data.Stochastic simulations are good for local calculation of the distribution of the studied variables, balancing the location b a last recording year analysed in this study).To facilitate comparison, the colour code for nitrate concentrations is the same for all five images.This sequence of images shows that the highest values occur in parts of the coastal plain of Valencia and Castellón, with one in the plateau of La Mancha; most of the interior areas are not monitored and thus the values between known data points are estimated with distant data.Also, to simplify interpretation of the zones with higher nitrate concentrations, binary images of the same hydrological years were computed (above or below the reference value of 50 mg/l), and are represented in of the data, the variograms and the contradictions between the true data.For this reason, to compute binary images of a variable above or below a threshold, it is preferable to use the set of simulated variograms instead of transforming an estimated image to a binary one.However, an average image obtained by a large number of simulations is almost equal to one obtained by kriging.Fig. 9 displays a simulated scenario of nitrate concentrations for the hydrological year 1970/71 overlaid with the measurement locations of the same year.This simulated image shows a possible scenario of nitrate distributions respecting the theoretical model of variograms fitted.
In the analysis of temporal evolution, Fig. 10 shows several examples of the results of the computed average images of nitrate concentrations for the hydrological years 1974/75, 1984/85, 1994/95, 2004/05 and 2011/12 (the

FINAL REMARKS
Groundwater contamination by nitrates from agricultural and agro-industrial activities in the Júcar River basin region reaches worrying levels in some areas, particularly the coastal plain areas and the plateau of La Mancha.Nitrate levels increased steadily from the early 1970s until approximately 2005, since when the values appear to be more stable.In the agricultural region of Albacete, in the continental plateau of La Mancha, where cereal crops and vineyards are concentrated, the situation is not so alarming, but nitrate levels are increasing.
The quality of water in the Júcar River basin is of particular importance because degradation reduces availability in this semi-arid region.Overall management and sustainability of surface and groundwater is essential to ensure a continuous supply for different uses, including urban areas, tourism, agriculture and industry.Interdisciplinary studies involving geology, hydrogeology, geochemistry and geostatistics are fundamental to management of water resources and to achieve desirable quality levels.

Fig 2 .
Fig 2. Digital terrain model showing the most important geomorphological structures.Fig 3. Groundwater bodies (GWB) of high permeability, and those of low permeability or local interest only, overlaid with sampling points with records of nitrate concentrations in the Júcar River Basin District.

Fig 6 .
Fig 6.Classification of groundwater bodies according to their dominant lithology.

Fig 9 .
Fig 9. Result of one simulated image of nitrate concentrations (mg/l) by sequential Gaussian simulation algorithm for the hydrological year 1970/71.
Fig. 11.Overall, both Figs. 10 and 11 show that the problems of contamination in mostly agricultural areas, such as the plateau of La Mancha and the coastal plain, are persisting and worsening.Nitrate concentrations appear to have increased steadily from the early 1970s and to have stabilized since 2005.
likely to create dissolution structures like karsts), and carbonate, detrital type of the Cretaceous period.Fig. 6 displays the groundwater bodies according to the classification of the rocks: (a) carbonate, crystalline type; (b) carbonate, detrital type; (c) detrital type; (d) low permeability or of local interest only and (e) not considered for this study.
detrital, carbonate or mixed, the last designation being applied if both detrital and carbonate materials occurred.Carbonate rocks were also classified in two subclasses: carbonate, crystalline type (mainly from the Jurassic period, more

Table 1 : Number of positive and negative correlations observed for the studied groundwater bodies (GWB) according to the type of predominant rock type and climate variable (rainfall or recharge from PRATICAL model) Correlation coeficient signal Seven groundwater bodies mainly of carbonated crystalline rock-type Rainfall Recharge calculated with PATRICAL model Same year One year of delay Two years of delay Same year
in groundwater bodies of carbonate lithologies and direct correlations in groundwater bodies of detrital lithologies.In carbonate rocks, the positive correlation coefficients can be explained by the prevalence of the detrital component over the crystalline component, or by the existence of a significant detrital component that overlies the carbonate formation.However, for detrital formations, the correlation between nitrate concentrations and rainfall or infiltration delayed by one or two years led to improvements in the correlation coefficients.As previously stated, and unlike with carbonate rocks, positive correlation coefficients predominate over negative ones.
Fig 7. Distribution of trends in nitrate concentrations in the studied groundwater bodies.correlations