Qualitative and Comparative Study of Different Methods of Interpolation for the Mapping of Groundwater Salinity: Case Study of Thermal Waters Used for Irrigation in Northeastern Algeria ()
1. Introduction
Northeastern Algeria, like many regions of the world, has experienced for the last two decades a decrease in the recharge of its hydro-systems. To cope with this, agriculture has gradually been increasing its abstraction from underground resources by drilling or collecting from springs. The geology of this region is mainly karstic in nature, and these waters rising from deep reservoirs have variable hydrochemical characteristics linked to the rocks they pass through and their residence times. However, over time the composition of water can have a detrimental effect on irrigated soils by progressively reducing their yields, particularly by sodization or alkalinization (Amroume, 2018).
This study proposes the use of a Geographic Information System (GIS) integrating a geostatic analysis model (Johnston et al., 2001) to map the salinity (Cemek et al., 2006), characterize the chemical quality of aquifer systems in North-Eastern Algeria, and provide information on the total mineralization of groundwater from karst springs (Feraga, 2017) that could be used for irrigation. In spatial statistics, the geographic positions of data are recorded so as to apply this spatial information in statistical modeling through cartographic representation (Allard, 2012). Kriging is a stochastic method of spatial interpolation that predicts a value (which can be the expression of a natural phenomenon) at unsampled sites, by an unbiased linear combination of minimal variance observations of the phenomenon in neighboring sites. In this study, a methodology for implementing kriging is proposed and illustrated. The kriging performances are then compared to the interpolation Inverse Distance Method (IDW) to estimate the risks to agricultural areas irrigated with mineralized groundwater. The spatialization method has already been compared in hydrology (Ahmadi & Sedghamiz, 2007) and agronomy (Amini & Khadi, 2002; Koussa, 2018) but never for deep groundwater at this regional scale in a karstic complex region, where tectonic faults control the location of springs with highly mineralized water. These hydrogeological and climatic contexts can be found in other parts of the world, including in North Africa, the Middle East, China, and also around the Mediterranean basin (Figure 1).
2. Presentation of the Study Area
In the North East of Algeria, several thermal springs with strong flow rates are located in a region of 100,000 km2 which extends over a width of 500 km from west to east, from Bejaia to El-Kala (border with Tunisia), and 200 km from North to South.
2.1. Physical Appearance
Successive orogenic phases throughout the different geological eras are very distinct and are responsible for the morphology of this region: 1) The coastal reliefs or Tellian area which include the Babor mountains to the west, the M’cid Aicha mountains in the center, and the Edough mountains to the east; 2) The high plains of Setif and Constantine, including vast flat and monotonous stretches, occasionally interspersed with chotts and sebkhas; 3) The southern massifs which include the mountains of Hodna, Bellezma, Aurès, Tebessa de Medjerda and Mellègue.
2.2. The Climate
The region has a Mediterranean climate. The coastal region and the highlands are characterized by constant rainfall in Northeastern Algeria (Nouaceur et al., 2013), and have experienced a slight annual decline mainly to the West since the mid-1970s, with the decade of the 1980s being the most deficient. Over a period of almost 80 years, the interannual minimum temperatures were around 16.5˚C on the coast and 14˚C in Constantine, in the center of the study area. The interannual monthly average temperatures vary between 7˚C and 28˚C over the whole area except in the highlands. However, the average temperatures have increased by more than 0.5˚C since the 1980s; with 1992-2002 being the hottest (about ±0.7˚C), especially in Northern Algeria.
The downward trend in precipitation and the increase in temperatures are primarily explained by the positive phase of the North Atlantic oscillation according to (Taibi, 2011). These changes have accentuated the use of groundwater for agricultural uses, in particular due to the decline of surface resources.
2.3. Geology
The geological context is marked by a geological structure set up by tectonic episodes (alpine orogeny) of Priabonian to Tortonian age. These events have created large faults almost with E-W and NE-SW directions and compressive structures such as folds and overthrusts. In recent periods, from the end of the Miocene to the present, extensive regimes induced movements of faults leading to karstification and the rise of deep water. The region is now rich in springs with thermal anomalies.
The area consists of limestone, marl, marly limestone, and sandstone that form the mountainous setting with some Triassic outcrops. There are seven main geological units: the Tellian unit, limestone chain, Neritic unit, South Constantinois, South Sétifian, parautochtone, and flysch (Figure 1).
The Quaternary period saw the creation of the Hamma Bouziane, Tamlouka de Bordj-Beïda, and Teleghma depressions corresponding to the set of faults that may still be active where several hot springs are located. These depressions are favorable to agricultural activities that require extensive irrigation.
In the north (along the coast) the area is characterized by small displacements frequently finishing in overlaps and outcrops of the crystalline bedrock which runs from the west to the east, often in direct contact with very large limestone masses (limestone unit chain).
To the south along the transversal sections, there are three successive types of structure organization from west to east, highlighted by diapirs from the Triassic (Vila, 1980) in the foreland Tellian groove: the setifian organization, Constantinian organization, and the Algerian-Tunisian organization. The Alpine overthrust is responsible for large overlaps and shear zones usually located at the base of the different units (Vila, 1980). It is also responsible for the anarchic appearance of detailed structures in the Tellian unit, and for the cutting of the folded Atlasic structures of the Autochthonous.
2.4. Hydrogeology
Northern Algeria was the subject of numerous hydrogeological studies (Energoprojekt, 2009) and localized geophysical prospections. These revealed that virtually all of the Triassic formations (gypsum and clay-salt) have very low permeability and cannot constitute real reservoirs; the springs which emerge nearby are highly mineralized with an evaporitic signature linked to dissolution, in particular of halite or gypsum.
However, the outcrops of the Jurassic, Cretaceous, and Tertiary periods which have a carbonated, marly, or sandstone nature, and are widespread in northern Algeria in particular around the Tell Atlas, constitute powerful reservoirs. The karst aquifers, which are the most common, are exploited by boreholes or drained by high flow sources. The geological structure of these reservoirs is characterized by folds and faults that favor the emergence of thermal water along NE-SW tectonic accidents (Feraga, 2017).
The late Mesozoic and early Paleogene eras, characterized by clay and marl sedimentation, produced levels with a very low water capacity.
Paleogene formations composed of sandstone conglomerate and limestone can constitute small reservoirs.
Quaternary sediments unconformably overlie the most ancient formations in the highlands. Thus, gravels, sands, silts, and gypsiferous clays present characteristics favoring important aquifer prospects that could be very useful for local agricultural uses.
Figure 1. Location map of Northeast Algeria (red box)—(from Nations Online).
This regional approach to the hydrogeology context shows that the karstic aquifers are the most important at this scale, with the most powerful springs.
2.5. Vegetation
The plant cover is composed of forest species (Sini perus, Pinus and Cedrus). This is the case for Djebels Abd-en-Nour Fortas and the Babors (Guerzouli, 2014). However, some calcareous areas or dolomitic mountains of Hodna and Constantine, where rainfall is of the order of 300 mm to 400 mm/year are comparable to real deserts due to the efficiency of infiltration.
2.6. Agriculture
In the north of Algeria, the agricultural area is characterized by the presence of two mountain ranges, the Tellian Atlas and the Saharan Atlas, stretching from the western border to the eastern border. Between these two chains, the high plains of Setif and Constantine are spreading. The majority of plains suitable for agricultural activity are marked by aridity or semi-aridity, and the majority of wetlands are mountainous. Their conjunction exists only in very limited areas, such as in Metidja or the plains of Annaba-Skikda in the east, where rainfall is greater than 600 mm and the slopes are less than 3%; these areas cover only 500 000 ha. Only 30% of irrigated surfaces receive more than 400 mm of water a year and the arid and semi-arid ground represents 85% of the total area (Hachaichi, 2018; Bessaoud, 2019). In addition, the mediocrity of the hydrographic network does not favor the development of an intensive agricultural area. In these conditions, springs with good permanent flow rates would be an evident solution for irrigation.
3. Materials and Methods
All data processing was carried out on ArcGIS version 10.3 with the Spatial Analyses extensions, based on the results of hydrochemical analyses carried out on 51 samples of karstic spring water.
3.1. Inverse Distance Weighted Interpolation IDW
Inverse distance weighted (IDW) is a local method of interpolation. In this method, the value estimated by a point in the study area is determined using the weighted average of the point values closest to the point considered. At the start, we measured the distance between the considered point and the known points in the vicinity. Following this, the sought point is calculated using the average of the values of the surrounding points. Thus, the more the point to be interpolated is close to a point of known value, the more the interpolated value will be close to the known value (Leroux, 2007):
where
P is a real positive number, called the power parameter. Here, the weights of the neighboring points decrease as the distance increases. The largest values of P give more influence to the nearest values of the interpolated point.
3.2. Ordinary Kriging
The hypothesis of simple kriging that the expectation of the random function Z(·) is known is rarely verified, and this method has therefore been generalized for cases where the expectation is unknown but locally constant, i.e. in the kriging neighborhood. Ordinary kriging (Matheron, 1970) is the kriging technique most frequently used following Gratton (2002). This type of kriging requires no stationarity assumption of order two. Thus, it will be developed here under the more general hypothesis of intrinsic stationarity. The basic model of this method is as follows (Baillargeon, 2005):
With μ unknown and almost constant and δ(·) a stationary random function intrinsic with null expectation and known dependency structure. The almost constant character of μ means that expectation is not constrained to remain the same all over the field D. On the other hand, it must remain constant within each kriging neighborhood (Arnaud & Emery, 2000). Thus, a forecast at point S0 is based on the following vector model:
with
Here is step by step the resolution of the equations of the ordinary kriging, based on (Cressie, 1993).
Linearity constraint
The forecast
must be a linear combination of the random variables
to
. It is therefore expressed in the form:
No-bias constraint
As with simple kriging, forecasting must be unbiased. It is therefore necessary that:
From the previous step, the constraint that
must be respected, thus, “a” must be set to zero, and the structure of the forecast therefore simplifies to:
with
Optimality constraint
The last step consists of finding the weights λi which minimize the variance of the forecast error. To do this, it is necessary to express this variance as a function of Г and
, as these statistics are known.
3.3. Sampling and Analysis
The hydro-geo-chemical approach undertaken is based on two sampling campaigns which were carried out within the framework of this study (Feraga, 2017); the first was carried out in the period of low water for all thermal springs in the study area, during September and October 2014, and the other was carried out in the period of high water in May, 2015. In total, 51 points were sampled (H. Azzaba was sampled in three wells very closed) and the results of the chemical analysis can be found in Table 1.
The analyzed parameters were:
Cations Ca2+, Mg2+, Na+, K+.
Anions Cl−,
,
,
.
Some trace elements.
Temperature, Electrical conductivity (EC), and pH.
EC is expressed at 25˚C in dS/m (equivalent to 1000 μS/cm).
3.4. Water Quality Assessment
Residual Sodium Carbonate Index RSC
To predict the evolution of the chemical composition of a solution, the concept of residual alkalinity was used (Barbiero & Valles, 1992). Carbonate and bicarbonate ions combined with calcium or magnesium will precipitate as calcium carbonate (
) or magnesium carbonate (
) under drought conditions. When the concentration of Ca2+ and Mg2+ decreases, the sodium content and the index RSC (Residual Sodium Carbonate) become comparatively greater. This will cause an alkalizing effect and a pH increase.
Therefore, when a water test indicates a high pH, it might be a sign of a high content of carbonate and bicarbonate ions. The RSC is calculated using the following equation:
· If RSC becomes positive (precipitation of calcite/sepiolite) following the addition of the sulfates relating to the precipitation of the gypsum, then it is referred to as a neutral saline route with a predominantly sulfate base.
· If RSC becomes negative even with the addition of sulfates, this is called a neutral salt pathway with chlorine dominance (Marlet & Job, 2006).
Irrigation water salt index SAR
A high quantity of sodium ions in water affects the permeability of soils and causes infiltration problems. This is because the sodium in exchangeable form in the soils replaces the calcium and magnesium absorbed in clay soils and causes the dispersion of particles in the soils (that is to say if the calcium and magnesium are the predominant cations absorbed on the soil exchange complex, the soil tends to be easily cultivated and has a permeable and granular structure).
This ion dispersion results in an alteration of soil aggregates. The soil then becomes hard and compact (when dry) thus reducing the infiltration rates of water and air, and hence affecting its structure.
This problem is also connected with several factors such as the salinity rate and the type of soil. For example, sandy soils do not damage so easily in comparison to heavier soils when they are irrigated with water with high SAR. The high levels of sodium become a problem when the infiltration rate is so reduced that the harvest does not receive enough available water, or when the hydraulic conductivity of the soil profile is too low to provide adequate drainage. Other problems for harvests caused by excessive Na+ (Bradaï et al., 2008) are the formation of crust beds, temporary saturation at the surface, high pH, and increases in the presence of grass diseases, soil erosion, lack of oxygen, and unsatisfactory nutrient availability. Recycled wastewater may be a source of excess Na+ in the soil compared to other cations (Ca2+, Mg2+, K+) and should therefore be appropriately controlled. The SAR can be calculated using the following equation:
SAR < 3 meq/L: no restriction on the use of water.
3 < SAR < 9 meq/L: between 3 and 6, special attention must be paid to sensitive crops and between 6 and 8; gypsum should be used, for non-sensitive crops.
SAR > 9 meq/L: severe damage.
Statistical analysis of data
Statistical processing allows the description of the data distribution and the statistical characterization of the studied variation (distribution law) and dispersion (coefficient of variation, variance). Figure 2 shows the location of the data (Koussa, 2018).
4. Results and Discussion
4.1. Water Quality Assessment
The sampled groundwater has been analyzed and presented in a Piper diagram (Figure 3). This makes it possible to characterize a body of water according to its geochemical facies, as well as to highlight groups of hydrochemical facies, by overcoming the differences in total mineralization between the waters.
Four water groups have been distinguished (Feraga et al., 2016):
1) 38% of the water is sodium hyperchlorinated (Cl−-Na+) and 2% sodium sulfated (
).
2) 32% of the water has a sulfated calcium facies (
) and 10% chloride calcium facies (Cl−-Ca2+).
3) 10% of the water has a calcium bicarbonate facies (
).
4) 10% of the water has a sodium bicarbonate facies (
).
Figure 2. Location of sampled springs and electrical conductivity of water.
Figure 3. Piper diagram of thermal waters in north-eastern Algeria.
4.2. Residual Alkalinity of Water
The evaluation of the quality of irrigation water (Durand, 1982) based solely on SAR very often minimizes the risk of sodization and alkalinization of the water when the chemical facies is not chlorinated. This is why we have considered the residual alkalinity as another approach to assessing the quality of the thermal waters of Northeastern Algeria (Table 2).
The use of each of the three types of water mentioned in irrigation presents a risk of degradation to varying degrees (Marlet & Job, 2006).
The RSC1 > 0 types of water present a major risk of alkalization and degradation of the physical properties of soils by sodization (Delbari & Amari, 2014), well explained by SAR with a value higher than 9 meq/L.
The RSC2 > 0 types of water mainly risk degrading the physical properties of soil following rapid sodization that SAR only imperfectly characterizes (3 < SAR < 9 meq/L).
The RSC3 > 0 types of water have no restriction (SAR < 3 meq/L).
The RSC4 < 0 types of water represent severe damage (SAR > 9 meq/L).
The RSC5 < 0 types of water exhibit sodization, and special attention must be paid to sensitive crops (3 < SAR < 9 meq/L).
For RSC6 < 0 types of water, the soil degradation risk is low since the sodization is progressive and salinity is then high enough to ensure the structural stability of soils (SAR < 3 meq/L).
These results indicate that the majority of the sampled groundwaters present a risk for soils irrigated with these waters.
4.3. Descriptive Statistics
The electrical conductivity (EC) of water is a function of its total mineralization or its salinity.
Over the 51 samples, the mean of the calculated EC is equal to 3.18 dS/m (Table 3). This value greatly exceeds the limit of 1 dS/m (Boufekane & Saighi, 2016) beyond which mineralization becomes high, and the tolerable limit of water intended for consumption is reached (Rodier, 1996). This is also the case for irrigation water, for which the limit is 2.25 dS/m, after which salinity becomes important and can harm plants and soil (USSL Staff, 1954).
Table 2. Percentage of each type of sampled water by RSC class.
Table 3. Statistical parameters of electrical conductivity of water.
The number of samples with an EC greater than 2.25 dS/m equals 21 samples, or 41.2% of all 51 samples. The class below 1 dS/m contains only 5 samples, against 6 samples with salinity above 5 dS/m.
4.4. Inverse Distance Weighted Interpolation IDW
Four threshold values for irrigation water were chosen: EC = 1 - 2.25 - 3 and 5 dS/m. 1 dS/m is considered in the literature to be an upper limit for moderately mineralized water, 2.25 dS/m is considered in the literature as a warning limit for irrigation water, 3 dS/m is considered as very high and the use should be restricted, and above 5 dS/m the water is considered as dangerous for plants.
We note that the interpolation methods (IDW, ordinary kriging) show similar results to each other, with irregular spatial distribution of EC in different areas (Figure 4, Figure 5).
The shape of the obtained maps reflects the fact that the methods used greatly depend on the location of the sampled points and are based solely on the neighborhood. We notice that the electrical conductivity maps show a certain spatial uniformity of the different levels of distribution of EC in the study area. It is interesting to notice that kriging shows structuring NE-SW lines that are coherent with the direction of the major faults on which sources are positioned.
The values of the electrical conductivity were adjusted to the experimental variogram in an exponential model with a nugget effect of 0 dS/m, a bearing of 41.98 dS/m2 and a range of 10,388 m (Figure 6).
The high value of the plateau also expresses a high spatial heterogeneity. The existence of a range justifies the use of ordinary kriging (Table 4).
4.5. Comparison between IDW and Ordinary Kriging
In this part, we have performed a comparative analysis between the two different methods (determinist and geostatistical) to choose the representative map of the distribution of conductivity in the spring waters of North-Eastern Algeria.
The adjustment of the point clouds measured between EC and residues estimated by IDW and ordinary kriging gives a straight line with positive slope, despite significant dispersion (low R2). This shows that the tendency is to overestimate rather than underestimate (Figure 7, Figure 8).
The statistics of the mean of the estimation errors and the standard deviation of the estimation errors of different methods are summarized in Table 5.
The results of Table 5 show that ordinary kriging is more accurate and less biased because the average of the errors is near to zero and the standard deviation of estimation errors is smaller.
Furthermore, kriging seems to reproduce better the organization of the groundwater resources since NE-SW alignments are visible. This direction corresponds to the large tectonic faults that allow deep water (often thermal) to flow upwards, giving rise to springs.
Table 4. Characteristics of the experimental variogram used for the ordinary kriging method.
Table 5. Statistics of estimation error for the different methods.
Figure 4. Distribution map of the electrical conductivity EC (dS/m) by IDW.
Figure 5. Distribution map of the electrical conductivity EC (dS/m) using ordinary kriging.
Figure 6. Model of semi-variogram adjusted to data and results of ordinary kriging for conductivity (Y axis: variance).
Figure 7. Fit of point clouds between measured ECs and residuals (IDW).
Figure 8. Fit of point clouds between measured EC and residuals (ordinary kriging).
5. Conclusion
Based on the sampling of 51 deep groundwaters coming from karstic springs, it was possible to compare two methods of mapping the electrical conductivity of water. This physicochemical parameter is considered to be representative of the salinity, which can alter the quality of the irrigated soils from the catchment of these waters and give the hydrochemical characteristics of the waters (RSC and SAR values). The kriging method appears to have some advantages over the inverse distance method (IDW).
Mapping of the electrical conductivity of groundwaters in northeastern Algeria has shown the advanced state of the water’s vulnerability to salinity. The estimation made by ordinary kriging shows that only 9.3% of the total area of the zone contains groundwater with a salinity of less than 2.25 dS/m, while 57.3% of this area is greater than 3 dS/m.
Finally, this study proposes an original approach to the consideration of the use of thermal karstic waters for irrigation purposes. These waters are very special due to their origin, which explains their composition. However, they are used because no other resource is available. From the point of view of groundwater management for irrigation, the long-term use of thermal waters poses a potential risk for soil salinization, and thus for crop productivity in Northeastern Algeria. This risk is accentuated by an arid and semi-arid climate in this irrigated zone, which increases the danger of salt concentration in the root zones of the soils. From an agronomic standpoint, it would be possible to carry out a control of the chemical quality of the soils, in particular their salinity or their SAR index, in order to correct it if necessary, for instance by adding calcium sulfate.
However, this approach could be further refined by finely analyzing the relationship between water electrical conductivity (EC) and the water quality indices RSC and SAR for several hydrochemical profiles. This could reinforce the robustness of using the EC of water as a proxy for the quality of water used for irrigation. In regards to mapping, a kriging-oriented approach according to geological structures could be relevant when it has been demonstrated that these structures play a role in the position of springs (NE-SW normal faults in the case of Northeastern Algeria).
This approach could also be used in many other arid or semi-arid countries where climate change is inducing a reinforcement of catchments based on groundwater from the deep origin in karstic reservoirs, providing salinized waters in springs or wells. A world map of sensitive areas could be drawn up by comparing tectonized karstic areas (favorable to the presence of thermal springs) and agricultural areas in arid, semi-arid, or even Mediterranean climates. With this kind of perspective, this study will provide an efficient stepping stone for the orientation of agricultural policies according to the quantity and quality of groundwater resources in many parts of the world.