GIS Applied to the Hydrogeologic Characterization – Examples for Mancha Oriental Aquifer (SE Spain)

1. Introduction

The population on planet Earth, according to FAO forecasts,will increase from 6 billion to 8.1 billion inhabitants in 2030 andwill coincide with an increase in water demandsto meethuman needs. Fresh water has ceased to be an inexhaustible resource to become a rather limited and scarce one.

Earth's hydrosphere has an approximate volume of 1.38x10¹⁰ km³ of water, which has remained virtually constant since its formation over3 billion years ago. This volume of water is distributed into four groups:

1. The vast majority is in the oceans, at 97.6% of the total (1,350x10⁶ km³),

2. In second place is solid water, in glaciers, at 1.9% (26x10⁶ km³),

3. Third is groundwater, with 0.5% of the total, which is 7x10⁶ km³, and

4. The remainder of water on Earth (0.03% of total) is divided among lakes (0.017), soil (0.01%), the atmosphere (0.001%), Biosphere (0.0005%) and rivers (0.0001%).

Ocean water is salt water and the glaciers are difficult to utilize because they are located far from major populated areas. Therefore, we find that groundwater is the largest volume of freshwater available to man. The volume of groundwater is 4,000 times greater than that of rivers and 30 times higher than the rest of liquid water that is on the surface of the continents. In addition, groundwater has characteristics that make it especially attractive to combating the processes of drought and desertification. Unlike surface water, groundwater does not evaporate, there are no major seasonal variations and flow is very slow, so itis difficult to contaminate (Castaño et al., 2008).

Groundwater is held in a naturally occurring reservoir called an aquifer, a geologic formation capable of storing, receiving and transmitting water so that man can easily take advantage of economically significant quantities to meet needs. The water is contained in any geological formation (ie: river gravel, karstified limestone, porous sandstone and so on).

Like all scarce resources, groundwater management must be approached from a dual approach (Knowledge and Sustainability):

Knowledge: There must be a sound understanding of the hydrogeological aquifer system to be managed. This should include a detailed analysis of the hydraulic aspects (geology, hydraulic parameters and groundwater flow) and should be contrasted with the hydrochemical aspects of water containing (origin of the substances dissolved in groundwater and hydro-chemical changes due to movement through groundwater flow).

Sustainability: Groundwater pumping of a region (an aquifer) shouldn’t exceed the water received, that is, the available resources. If users pump a volume of groundwater for short-term needs (ie: drought conditions), beyond the resources available, they use the aquifer reserves. Then, the aquifer must be given time to recover (either by saving water or allowing recharge to increase during periods of more rainfall). Otherwiseresources suffer overexploitation, putting the aquifer at risk of becoming depleted. It is obvious that the water volumes involved should be determined as precisely as possible.

In managing groundwater resources, Geographical Information System (GIS) are tools capable of storing and managing spatial hydrogeological data by spatial referencing in digital formats. The correlation of all data with location is the key feature of GIS, which provides the ability to analyze and model hydrologic processes and produce results in maps and in digital formats. Thus, GIS can be considered a support systemin decision making and an ideal tool for monitoring certain hydrogeological processes with socio-economic impacts (Goodchild et al., 1996).

Figure 1 shows a diagram of a GIS aquifer system modeling tool (Case study of the Mancha Oriental System). This scheme is integrated into a) a block of hydrogeological data maps of the study area, which supplies data on groundwater from urban and industrial and general hydrological information (surface and groundwater hydrology) necessary to carry out the integration and interpretation of some results, and b) a block of data from remote sensing imagery. Remote sensing allows for classification of crops and their relationship to water supply for irrigation, mapping of wells, and the assessment of recharge by precipitation of rainfall. All this information is transferred to software that simulates groundwater flow in the Mancha Oriental aquifer using intersection tools.

The main goal is to show the methods some GIS applications have in hydrogeological studies. This chapter is divided into two sections which describe some examples of hydrogeological characterization, and secondly, a method for calculating groundwater abstraction. To demonstrate these applications one of the largest aquifers in southern Europe (in terms of area), the Mancha Oriental System has been chosen.

Figure 1.

Basic diagram of coupling remote sensing and GIS techniques with the groundwater flow model of the Mancha Oriental System.

2. GIS & hydrogeological characterization

In general, among the Earth Sciences and particularly hydrogeology, sources of data tend to be from points (wells, points of water, lithological columns, etc.) defined by a geographic location (UTM or geographic coordinates) and attributes (topographic top or bottom of a hydrogeologic entity, groundwater level, hydraulic parameters, concentration of a chemical compound). This type of data, usually measured in the field, must be spatially distributed in a continuous manner such that a value is given for any point within the space. To achieve this, interpolation or spatial estimation is used. This method derives an interpolation function that provides estimates for a point in space based on the points measured. GIS tools have incorporated algorithms which perform these operations with discrete entities (vector) and generate spatially continuous entities (raster, line models, etc.). In addition to expanding and the geodatabase and adding values, these techniques create a foundation for spatial modeling (Peña Llopis, 2006).

The most commonly used spatially continuous entities are raster maps, which are characterized by a two-dimensional numerical matrix or digital image. Each element of the matrix, called a picture element or pixel, has an attribute assigned to it in the database. The only requirements are for maps to have attribute values referenced to the same coordinate system and the same number and arrangement of pixels to perform algebra operations with them (ie: isopaches: difference between top and the bottomraster maps of the geologic formation; calculation of storage volumes: difference between raster maps and contour lines or groundwater for different dates, multiplied by the storage coefficient, etc.).

2.1.Theoretical foundations

Using a spatial domain and a series of points (which we will refer to as points of observation) Pi, where i= 1, 2, ….n, which have a series of coordinates xi where variable Z has been measured in Pi points; Zi (observed values), the interpolation or spatial estimation aims to find the value of Z (estimated values) at any point in the known space. An interpolation function must be obtained:

Z(x)=θ(x,xi,Zi) =f(x)

The interpolation function should have certain characteristics: a) accuracy: the estimated value in points of measurement should coincide with the measured value, b) spatial continuity, c) ability to be derived: the interpolation should be “smooth” and d) it should be stable with respect to the location of the variable as well as its value such that small variations in data do not provoke large variation in the interpolation. As a function of these characteristics, especially the last condition, there is no universal interpolator and there is always another interpolation method which can be applied (Samper & Carrera, 1996). Most GIS software presents two interpolation methods: Deterministic and Stochastic.

2.1.1.Deterministic methods

This type of method is characterized by associating a mathematical function, such as an interpolation function, to the measured or observed values, in which these points are considered without error. Following the nomenclature followed until now, this mathematical function could be written in the following manner:

$$Z(x)=f´(x)={\displaystyle \sum _{i=1}^n{c}_{i}f{(x)}_i}$$E1

where for each x a Z(x) value is measured through a function f(x), which is defined by the sumof all “n” points of observation of a product between a base function f(x)i and coefficients, Ci. For example, in a simple exact interpolation the observed or measured values (Zi) coincide with the Cs values, multiplied by a weighting factor given by the function f(x)i. The deterministic interpolation functions differ from one another in the means of evaluating f(x)i and Ci.

There are various deterministic interpolation techniques. The most commonly used methods are presented here (ESRI, 1997) (Samper & Carrera, 1996):

Nearest. neighbor (Thyessen polygons, Polygons of influence).

This method assigns the value of each measured or observed point to each pixel or node of the interpolated area. For each point of observation the Euclidean distance is calculated for all other points and each is given the closest value. The result is a map of polygons with an interpolated value (Fig. 2A). This method is often used for regular grids and/or dense observed data, or to find areas of influence.

Interpolations. based on weighting functions.

The estimation or interpolation in this type of method is performed by a weighted average of the observed values. At each point of observation a weight is assigned. The selection criterion is that the weighting function is exclusively dependent upon the distance (d). The weight will decrease with increasing distance between points. The most common strategy for generating this criterion is the Inverse of Distance raised to some exponent (a).

$$f(x)=\frac{1}{d^a}$$E2

This exponent shows the “speed” with which the weight of a point of observation decreases with distance from the point of estimation. At times the number of points of influence is restricted, or a radius or maximum distance is assigned for considering points of observation. This interpolation method is exact and is commonly applied, with the only disadvantage being creating the feared “bulls-eyes” (Fig. 2B).

Polynomial. interpolation.

In this method the interpolating function is a polynomial function which varies in its exponential order. The choices for polynomial are: a) through exact fit and b) fit by least mean square. The first method aims to resolve the system of equations defined by the n points of observation. If there are many points of observation, the fit of higher order polynomials can become unviable, giving unrealistic interpolations with exaggerated variation among the values (Fig. 2C). In fact, by default these methods limit the polynomial to third order and only use the number of points in a nearby group. One special case of polynomial interpolation is linear interpolation, wherein the interpolation function is a first order polynomial which directly depends on the position of the observed values. It is an exact method and does not take into account the spatial distribution of the variable, with the result of soft surfaces. It is an easy method and is often used, above all in cases when not a lot of data is available and the aim is to study the spatial variation of a certain variable. In general, this interpolation method is not used for spatial estimates on realistic structures (topography, groundwater levels, etc.) but rather to determine the tendency of data (Fig. 2C).

Spline functions.

Within polynomial interpolation, this general method generates a different series of expressions for each subdomain into which the whole interpolation space has been divided, wherein continuity requisites are imposed, especially in the contours common to more than one subdomain. The results of this interpolation tend to be surfaces with small changes in levels (Fig. 2D).

2.1.2. Stochastic methods

This methodology is based on the premise that the variable to be interpolated is a random function associated with probabilistic distribution laws. This type of method gives a

Figure 2.

Examples of the different deterministic interpolation methods using GIS tools. A) Polygons of Influence (Thyessen polygons, Nearest neighbor), B) Interpolations based on weighting functions (in this case the function is the inverse distance squared (IDW); C) polynomial interpolation, in this case third order; D ) Spline functions. The same data are used for each case and the results are shown as continuous identical spatial entities,with the same size and pixel size.

measure of error of the interpolation based on the data. There are two classes of stochastic interpolation: a) non-parametric, which are not exact because the errors are assumed to be independent and b) parametric, wherein the interpolation function depends on certain parameters calculated as a function of the observed data (IDW or Krigging) (Samper & Carrera, 1996). The most common method, Krigging, which is available in most GIS software packages, is explained below.

Krigging was created under a new discipline, geostatistics, as a resultof problems presented by deterministic interpolation in Earth sciences due to the uncertainty and variability of data (Cassiraga, 1999). The starting hypothesis of geostatistics is that the data of study has a correlation spatial structure, as the realization of an infinite amount of possible realizations. For this reason, geostatistics is called the science of regionalized variables. For spatial estimation using Krigging, the steps below should be followed, among others, and the variable to be interpolated should meet the criteria of normality and stationarity (Johnston et al., 2001).

The first step is structural analysis, with the objective of estimating the semivariogram. This relates the Euclidean distance among the points of observation with the variability of the measured values (Samper & Carrera, 1996). First, the variable should be defined as stationary, if there is a tendency among the data, etc. The function of the semivariogram (estimator of spatial variability) is expressed as:

$$\gamma (x)=\frac{1}{2N}{{\displaystyle \sum _{i=1}^n\left[Z\left(x_i+h\right)-Z{\left(x_i\right)}^{}\right]}}^2$$E3

where:

Z(xi): experimental data,

h: distance between points of observation (Variogram step),

N: number of pair separated by vector h found in a group of data,

Xi, xi+h: experimental points in an n-dimensional space.

At first, from the observed data an experimental semivariogram will appear. A theoretical function with similar behaviourcan be fit to this in order to calculate a weighting matrix for each point, and statistical error affecting the interpolation can be calculated. The semivariogram is composed of a series of elements (Fig. 3A):

Range: the distance from which the spatial correlation is practically null (Area of influence),

Sill: value that the semivarogram takes in the Range,

Nugget: value of the semivariogram when it intersects at the coordinate axis.

Figure 3.

A) Elements of a spherical semivariogram. B) Experimental and theoretical semivariogram sets.

The experimental variogram cannot be used for the geostatistical application. It must be fit to a theoretic model (Fig. 3B). There are different technical variogram models available, with the most popular being the stationary or spherical semivariogram. Once the theoretical semivariogram has been chosen, the Krigging technique performs the spatial estimation of the data. There are diverse Krigging techniques as a function of diverse methodological hypotheses:

Simple: Hypothesis of stationary variable with a known mean and covariance,

Ordinary: Hypothesis of stationary variable with anunknown mean and known covariance,

In an environment (by blocks): Quasistationary hypothesis,

Residual: Non-stationary hypothesis with a known drift, from which residuals are derived and ordinary Krigging can be performed,

Universal: Non-stationary hypothesis and polynomial form with a drift set a priori.

2.2. MOS case study

The Mancha Oriental System (MOS) is located in the SE of the Iberian Peninsula and is one of the largest aquifers in Spain (7,260 km²) (Fig. 4). The area has a semiarid Mediterranean climate. Average rainfall is 350 mm/year and mean annual temperature is 13-15 C; the continental nature of the climate is clearfrom the extreme temperatures that occur.

The area is characterized as a high plain (700 masl mid-altitude) surrounded by gentle relief, interrupted only by a valley which was carved by the action of the JúcarRiver. From a hydrogeologic perspective, the MOS is formed by the superposition of three limestone aquifer hydrogeologic units (UHs): UH2: Tertiary, UH3: Upper Cretaceous and UH7: Middle Jurassic. These HUs are separated by aquitards/aquifuges that comprise UH1 (upper and lower), UH4, UH5 and UH6 (Sanz et al., 2009). The impermeable base and southwest boundary of the area of study is composed of marl, clay and gypsum from the Lower Jurassic, belonging to HU8 (Fig. 4).

Over the last 30 years the progressive transformation of approximately 100,000 ha from dry to irrigated farmland has translated into an acceleration of socioeconomic development due to widespread use of groundwater resources. Groundwater abstractions in the MOS exceed 400 Mm³/yr, of which 98% is used for irrigated agriculture and the rest to supply a population of 275,000 Inhabitants (Estrela et al., 2004). Groundwater pumping is not sustainable with the amount of available resources, estimated at 320 Mm³/yr by the Júcar Water Authority. Therefore, two major impacts are occurring: (a) the quantity of available groundwater is descending, noted as a continuous decline in the regional water level and a decrease in aquifer discharge to the JúcarRiver; and (b) the quality is also affected, as researchers have found a significant increase in nitrate concentrations in groundwater (Moratalla et al., 2009).

In this context, the MOS is an ideal case study for testing and validating the usefulness of GIS Techniques for understanding the aquifer system and planning for sustainable management.Following is a description of the interpolation methods applied to these variables: a) The elevations of the top and bottom of the aquifer units, b) Hydraulic parameters, c) Groundwater level data d) Groundwater chemistry. The approach is to explore the variable data with histograms and spatial trend analysis in orderto understand the behaviourof the variable in space as well as to establishwhether the data are consistent or anomalous. After analyzing the data, a variable is selected to be interpolated. The type of interpolation to apply for raster maps or continuous spatial entities is chosen.

Figure 4.

Simplified hydrogeological map of the MOS. Taken from Sanz et al. (2009). JRBJucarRiver Basin.Graphical output from GIS software.

2.2.1.Hydrostratigraphic framework

Any attempt at making a coherenthydrogeological model should be approached by first understanding, with a certain amount of precision, the geometric configuration. In addition to information on the surface geology, lithologic columns should be analyzed and gathered from water sampling points, and tested materials should be classified within the defined hydrogeologic units (Murray & Hudson, 2002). Using this information, layers of geographically located points (X, Y coordinates) as well as the topographic height of the superior (top) and inferior (bottom) limits were made into attributes of each hydrogeologic unit that behaves as an aquifer. Using geostatistical interpolation models developed on theoretic and applied foundations, GIS software (i.e. ArcMap ® 9.3) was used to determine the continuous geographic entities, ie: raster maps of the surfaces corresponding to the top and bottom of each hydrogeologic aquifer unit. The result is the three-dimensional structure of the hydrogeologic system (Fig. 5A). These 3D geologic models (Fig. 5), constructed using GIS tools, became the foundation for the numeric simulation models in later steps.

2.2.2.Hydrodynamic characterization

Transmissivity, permeability and storage coefficient are hydraulic parameters that must be quantified for an aquifer because they are neededto estimate the progression of groundwater levels, groundwater flow through a section of the aquifer, contaminant transport time, the degree of aquifer homogeneity and the numeric parameterization of the groundwater flow models(Mace, 2000). Generally, estimating these parameters requires pumping tests in specific points. These provide specific geographic entities defined by their coordinates, and attribute values for the hydraulic parameters in the well. It is also useful to have previous knowledge on the spatial behavior of the variable and establish a relationship for the interpolation model (ie: the type of distribution function of the variable). In the case of the Mancha Oriental System, to determine the spatial distribution of any of the parameters mentioned, these logarithms have been used because the variable tends to have a log-normal distribution. In this case, the value estimated by the Krigging method is the absolute optimum and the semivariogram better represents the structure of the spatial variability (Samper & Carrera, 1996). Once the structure of the spatial variability of log-T was studied, ordinary Krigging type interpolation models were applied (Fig. 5B and 3A).

2.2.3.Characterization of ground water flow

As is the case with aquifer hydraulic parameters, data on the height of the groundwater levels are also point data. The attribute of the groundwater level in this point geographic entity is compiled in the inventory of water points by subtracting the topographic height of the point from the depth of the water level in the well. These measurements should be performed for a specific date and in the least amount of time possible.

Figure 5.

Examples of hydrogeologic variable interpolation. a) The elevations on the top of HU7Mid-Jurassic aquifer, b) Hydraulic parameters (Transmissivity of HU2 Miocene Aquifer,) c) Groundwater level data for 1975 d) Groundwater chemistry (Nitrate rate-2008).

The raster maps obtained using the data on groundwater level height are called groundwater contour maps (isopiestic lines). These maps serve to determine how groundwater flow functions, where are the recharge and/or pumping areas in addition to indicating gradient calculations, flow and permeability (Fig. 5C). By crossing groundwater contour maps for different dates, contour descent maps can be obtained for that period. In addition, Variation in Saturated Thickness (VST) can be calculated between those dates.

2.2.4.Hydrochemical characterization

The chemical composition of groundwater is conditioned by a multitude of factors. Among those, the most important are: a) chemical composition and disposition of materials with which the water is in contact, b) time of contact with these materials, c) temperature, d) pressure, e) presence of gases, and f) level of water saturation in relation to distinct incorporated salts (Custodio & Llamas, 1983). Although the composition of groundwater is continually changing, the anthropogenic factors can significantly influence the composition. In fact, changes in land use are considered the most influential factor in groundwater pollution. Ions such as NO₃, SO₄, Na and Cl can come from agricultural fertilizers, livestock waste and waste from industry and urban centers. Nitrate is accepted as the most common contaminant in groundwater (Gulis et al., 2002; Jalali, 2009).

In Europe, the objective is for waterways to achieve “good” chemical and ecological status according to Directive 2006/118/EC of the European Parliament and Commission (DOCE, 2006). This directivedescribes the protection of groundwater from pollution and deterioration and the establishment of a pollution prevention and reduction plan by 2015. In addition, water bodies should be in good quantitative and qualitative status, especially in reference to nitrate content, which should not exceed 50 mg/l.

Thus, establishing the spatial distribution of NO₃ concentrations in groundwater within the aquifer is of vital importance. To accomplish this, point analyses of groundwater in wells and springs must be performed. Using advanced interpolation capacities provided by GIS tools, a complete geostatistical study can be performed to establish the most contaminated areas in terms of nitrate (Fig. 5D).

3. GIS & groundwater abstractions

Intensive use of groundwater for irrigation in arid and semiarid regions has often been the main driver of socioeconomic development over the past four decades (Shah, 2005). However, poor management of pumped volumes of water has led to negative consequences in terms of quality and quantity of available groundwater resources and associated ecosystems.

Controlling the groundwater withdrawals from a wide area of intensive irrigation is not easy. The largest volume of water used for agriculture has been extracted through tens of thousands of pumping-wells which generally have no measurement system and, in many cases, do not meet legal requirements or are unknown even in their location. Various methods of calculating groundwater abstractions have been known for years, but all of themare very expensive or inaccurate in their application to large areas (Brown et al., 2009).

In this scenario, the data provided by satellites (remote sensing) and the computerized processing of these geo-referenced data (GIS) represent a new approach to monitoring and quantifying groundwater abstractions, with the following characteristics: instantaneous observations are available over large areas, there are several images throughout the year, there is information not visible to the naked eye, data distributed in both space and time is available, the information is not conditioned by the legal or administrative characteristics of the pumping wells, and satellite image acquisition and processing is very low-cost compared to traditional methods (Castaño et al., 2009).

3.1. Theoretical foundations

The methodology for determining groundwater pumping for irrigation follows these steps: First, the irrigated crops are identified and classified by the multitemporal analysis of images obtained by multispectral sensors on satellite platforms, comparing the phenological evolution of the crops with the evolution of the Normalized Difference Vegetation Index (NDVI; González-Piqueras, 2006). Then, the area covered by crops is quantified by introducing the data into a Geographic Information System (GIS) and overlay them with the areas or required limits. Based on the surface area of each crop and the knowledge of their water requirements, the theoretical amount of water needed for those crops to reach the stage of development seen in the images is calculated. When agricultural practices are known, a correction factor is applied to translate the theoretical amount of water applied to each crop. Finally, all the information generated is integrated (spatially and temporally distributed) in a Geographic Information System (see Figure 1) and is used to establish relationships among all elements of the water balance (Brown et al., 2009).

3.1.1.Multitemporal analysis of satellite images and cross with vector cartography

The term “Remote Sensing” has different definitions, but the most commonly used is “a group of techniques that analyze data obtained by multispectral sensors located on airplanes, spatial platforms or satellites.” The sensors (on satellites) that observe the surface of the Earth are instruments that register the radiation from Earth and the atmosphere and transform it into a signal that can be managed in analog or digital format (Calera et al., 2006). The sensors do this by detecting the electromagnetic signal from the Earth and the atmosphere of a certain wavelength and converting them into an established physical magnitude. The energy values detected, quantified and coded from the sensors are usually in a two-dimensional number matrix or digital image (raster). Each element of the matrix, called a picture element or a pixel, has a digital value assigned to it (digital levels) which is usually registered in a byte or binary code (2⁸ values, from 0 to 255). These represent the energy associated with the wavelength to which the detector is sensitive.

According to Chuvievo (2002), each satellite scene can be used to extract four types of information, each with its respective resolution (Table 1):

1. Spatial, derived from the organization and presence of elements on the surface of Earth in three dimensions,

2. Spectral, dependent upon the observed and measured energy,

3. Temporal, associated with changes over time in a specific spatial location and

4. Radiometry, related with the conversion of voltage collected by the apparatus that receives the signal sent from quantifying entities and later on digital levels.

Landsat Satellite TM sensor				NOAA Satellite AVHRR sensor
Spatial resolution, pixel size	Spectral resolution	Temporal resolution	Radiometric resolution	Spatial resolution, pixel size	Spectral resolution	Temporal resolution	Radiometric resolution
30 x 30 meters	7 bands	15 days	8 bits	1x1 Km	5 bands	12 hours	10 bits

Table 1.

Type of resolution of two kinds of sensors, such as Landsat and NOAA satellite. (see http//:www.landsat.org and www.noaa.gov/).

The information referred from the sensor is treated digitally to obtain a geo-referenced representation of the land. Once the interactions of the atmosphere are removed from it, the radiation values received correspond exactly with those measured on the surface (see more detailed information in Chuvieco (2006) or Calera et al., (2006).

The source of radiant energy is solar radiation on the land surface after moving through the atmosphere. The radiation the sensor obtains is that which emerged from the land surface to the proper region of the spectrum when the emissions due to temperature are considered null. Therefore, the electromagnetic spectrum is the continuous succession of these frequency values (wavelengths). Conceptually it can be divided into bands in which electromagnetic radiation has a similar behaviour (Fig. 6).

Figure 6.

Main features of the electromagnetic spectrum.

Three basic elements can be distinguished as the components comprising all forms of the landscape on the Earth’s surface: soil, water and vegetation. The behaviour of these elements in different regions of the electromagnetic spectrum can be observed in Figure 7. The energy emitted (reflectivity) from the ground in the solar spectrum has a uniform response, showing a flat curve and ascending to greater wavelengths. It is important to know that bare soil can present different curves according to the chemical composition, humidity content, organic material content, etc. In the optical spectrum, water can be observed as a strong contrast between the reflectivity of the visible (5%) and the infrared, where water absorbs almost all the radiation in these wavelengths (Fig. 7). This effect is used to separate the water-soil limit. Similarly to soil, the characteristic reflectivity curve for water can vary as a function of factors such as depth, suspended materials, roughness, etc. (Calera et al., 2006).

Figure 7.

Reflectivity curve characteristic of soil, water and vegetation.

The morphology of the reflectivity curve against the wavelength that vegetation has is well defined (Fig. 7). It has low reflectivity (10%) with a maximum relative to the region of green, high reflectivity in the near infrared which is gradually reduced to the middle of the infrared spectrum. The strong contrast between the reflectivity of red and near infrared indicates that the higher the contrast is, the more vigorous the vegetation is, either due to greater land cover or greater photosynthetic activity (Calera et al., 2006). This spectral behaviouris the foundation for the development of certain indices with an objective to highlighting active vegetation from other components (soil, water, dry vegetation). From the reflectivity of each band (quantitative information distributed and geo-referenced in space) a relationship with the biophysical characteristics can be established (biomass, fraction of plant cover, etc.). This allows for quantitative, spatial-temporal monitoring of the processes on the Earth’s surface (Bastiaanssen et al., 2000; Calera et al., 2001; González-Piqueras, 2006). Nonetheless, reflecting the spatial and temporal variability of plant cover is complicated if different spectral bands with the reflectivity values are used. To unify this process, the Vegetation Indices have been developed, one of the most important being the NDVI (Rouse et al., 1973).

The NDVI is:

$$NDVI=\frac{(NIR-R)}{(NIR+R)}$$E4

where:

NDVI: Normalized Difference Vegetation Index,

NIR: Near Infrared reflectivity (spectrum range in micrometers),

R: reflectivity in red reflectivity.

GIS tools have the capacity to establish dynamic processes if they contain spatially referenced information which is repeated over time in addition to the ability to study spatial changes over the land surface.Mathematical operations can be used between the different sensor bands(digital images) to obtain quantitative information of each satellite scene obtained for a specific date.In this way, a temporal series is available for establishing the progression of a variable, for example NDVI.Multitemporal analysis stems from the availability of a time sequence of images, so these scenes must meet a set of requirements such as geometric coregistration (ability to superimpose images with the highest precision possible and radiometric normalization; Calera et al., 2005).

With this information and digital classification tools each pixel of the image from each date can be assigned a class defined through an automated process.There are two methods for classification: a) supervised and b) unsupervised (does not require intervention of an “interpreter”).The difference between the two is in the method of obtaining the spectral reference classes for assigning one to each pixel.

Supervised classification stems from a priori knowledge on specific land uses located in space, which are called training plots. These serve to establish spectral reference classes. There are several methods and a procedure for assigning a class to each pixel, but the most commonly used is an algorithm of maximum probability. Without getting into the details of this method, the algorithm is based on multivariate statistical analysis of components that identify each pixel in terms of their closest resemblance.

Other classification methods that could be used as alternatives or complimentary methods are decision tree (expert systems). These procedures are based on separating the pixel values of a layer into homogeneous groups and subgroups. Another method, called contextual filters, can also be applied. This not only considers the spectral characteristics of an individual pixel, but also considers neighboring pixels (Calera et al., 2006).

Once the classified map has been obtained, the spectral classes can be used to select the classes that are of interest from a hydrogeologic point of view.In our case study, this is crops irrigated with groundwater.Therefore, it is important to know the area of irrigated crops and their spatial distribution.One of the most commonly used GIS techniques for this is overlay vectorial and raster cartography, which is the only way to obtain this information for rasterized areas.There are two types of overlays that dependon the pixel value (ESRI, 1997). When the pixel has a real value (for example a precipitation map, groundwater level, NDVI, etc.) a statistic is calculated for areas by obtaining statistical values from the raster within the selected polygons (mean, minimum, maximum, etc.). The other case is when the attribute of each pixel is a discrete value defining a series of classes (i.e. raster map of the classified land uses). In this occasion (tabulate areas) the result is the surface of each class within the vectorial polygons selected (ESRI, 1997).Once the areas of each irrigated crop have been determined and the amount of irrigation water supply is known, it is possible to calculate the volume of water used to irrigate crops in the area on an annual basis (see summary in Castaño et al., (2009).

3.2. MOS case study

The use of groundwater resources of the MOS above its recharge capacity has led to several quantitative impacts: a steady decline in groundwater level, reduced aquifer discharge to the Júcar River and aquifer pollution. In fact, the quantitative analysis performed on the Júcar River Basin (Estrela et al., 2004) for the European Water Framework Directive (EU, 2000) clearly indicates that the environmental objectives set are not being reached at the present time and there is a certain risk of not meeting them by 2015.

In this situation, quite common in a semi-arid river basin, it is particularly important to precisely quantify the groundwater balance in order to determine aquifer sustainability. The information provided by the multispectral images becomes critical because these data sets are the only consistent and objective information on crops and can replace the data on agricultural statistics. In this regard, the MOS is an ideal case study for testing and validating the adequacy of remote sensing and GIS techniques for calculating groundwater abstractions in agricultural basins in semi-arid climates (Castaño et al., 2009).

Following is a description of several studies in the MOS to classify irrigated crops and quantify the ground water consumption required for ideal phenological development.

3.2.1.Calculation of groundwater withdrawals

The development of a method to calculate groundwater abstractions has been briefly described in the section on the theoretical foundations. In addition to knowing the method, one must have previous knowledge of the study area in order to choose the type of satellite image, for example crops and natural vegetation (phenologic development), soils, climate, relief, etc. In this study, due to the characteristics of the study area, the ideal images for thematic cartography are those from Landsat5 TM and Landsat7 ETM+ (Tables 1 and 2).

	Band 1	Band 2	Band 3	Band 4	Band 5	Band 7
Wavelength; (λ μm)	0.45	0.52	0.45	0.52	0.45	0.52
Region	Blue	Green	Red	Near Infrared	Middle Infrared	Far Infrared

Table 2.

Bands and wavelengths from the Landsat 5 TM and Landsat 7 ETM+.Information on band 6 from this sensor is omitted because it is not relevant in this study.

Using this information, the number of scenes required can be established as well as the bands to use for differentiating the crops of interest. For example, at least two images are necessary to establish a time series and to identify the non-irrigated crops, one in May or June (maturation process) and another in July (harvest). If more Landsat scenes on specific dates are included spring irrigated crops can be identified. Generally, a minimum of 16 images are used for performing the classification. If the temporal progression of the spectral response is considered a discriminating element for crops (phenologic development), the NDVI spectral band is the most useful. This is obtained by performing mathematical operations with the images for the same dates using bands 3 and 4 of the Landsat sensor (Table 2 and Figure 8).

Figure 8.

Example of the RGB composition for a maize plot and the temporal progression of NDVI during the crop growth cycle in the year 2005. (From Calera, et al., 2005).

The next step is to choose the classes to use in the classification as a function of those which can be differentiated using the spectral band in the images used such that they meet the objectives of the study (Spring Irrigation, Summer Irrigation, Spring-Summer Irrigation, Alfalfa, Bare soil, Dry farming crops, Shrubs, Forest).For classification, a study is required on the training plots called “true land.”

This information is used to perform the classification using the maximum probability algorithm, tree decision and contextualfilters. These classifications will be capable of discriminating between the sources of error that can be generated, e.g. in dry farmed crops because several images can have similar spectral responses. For example, in a rainy spring, cereals grown in irrigated or dry farming can be difficult to differentiate. Therefore, using other scenes with tree decision can help in the classification process. Contextual filters can be used to eliminate error on the plot boundaries or isolated pixels that belong to a different class than the rest of the plot. The result is a raster map with the classification of all irrigated crops (Fig. 9).

Figure 9.

Cropland irrigated with groundwater map, classified by remote sensing techniques for a2001 in the MOS from LANDSAT-TM scenes.

Once the maps have been classified, the spectral classes can be used to select the classes of interest from a hydrogeologic point of view. The area of irrigated crops must be determined (divided into spring, summer and spring-summer irrigated crops) as well as their spatial distribution (overlay tools). Information on the irrigation needs of the crops present in the MOS are provided weekly through the Irrigation Assessment Service (SAR) by the Agronomic Institute of Technology of the Province of Albacete (ITAP) using the method proposed by Allen et al., (1998). For each agricultural year, the institution groups the irrigation needs for each crop and publishes them in the annual monitoring reports (http://www.itap.es). The theoretical irrigation needs represent the minimum water consumption for sustaining the crops of interest. To determine the true water needs, correction coefficients must be applied to the theoretical irrigation volumes (Castaño et al., 2009). Field work to quantify the agricultural practices applied in the region should be done to perform this calculation. Thus, the irrigation efficiency can be used to calculate the correction coefficient that transforms the theoretic quantity of water necessary into the true values applied to each crop in the area. Generally the true amount of groundwater abstraction is higher than the theoretical irrigation needs.

Therefore, the calculation of water consumption for the different types of irrigated crops by applying the following equation:

$$Vr=A_i\times D$$E5

Where:

Vr = the annual volume of irrigation water for each type of crop (m³).

A = the area of each type of irrigated crop (ha).

D = the irrigation needs for each type of crop, applying the correction coefficient (m³/ha).

i = Hydrogeologic Domain, Municipality…In this way groundwater consumption could be calculated for the MOS or any part of it.

Estimating the amount of water required for irrigation is critically important in times of water shortage and especially in the current situation of increasing water demands with increasing populations. The use of GIS tools in this endeavor greatly increases the accuracy and efficiency of these types of study. This chapter is meant to be a summary of some methods used in the case study of the Mancha Oriental System, but they can be applied to other systems worldwide at risk of groundwater overexploitation or as a preventative measure to protect natural resources.

Acknowledgement

This study was funded by the Spanish Government under research grant CGL2008-06394-C02-02/BTE. Special thanks go to the Jucar Water Authority (CHJ) and Stakeholders (JCRMO) in the Mancha Oriental System for providing information. The findings presented therein belong to the authors.Special thanks go to Dr. A. Moratalla (UCLM) for providing some of the chemistry data and toDr. A. Calera and Dr. M. Belmonte (UCLM), for providing the multitemporalclassification of the crops in the study area. We would also like to thank S.A.Kroll for improvingthe English text.