ABSTRACT

This work focuses on the modeling and evaluation of water resources in complex aquifer systems and the use of scarce data. The modeling work is developed based on the GLUE (Generalized Likelihood Uncertainty Estimation) method. This method is still little used in hydrogeology, although its applications in other disciplines such as hydrology proved quite efficient. The study site, located in the Republic of Djibouti (Horn of Africa), is represented by the volcanic massif of Goda. The hydraulic properties of this massif are highly heterogeneous since they are associated with fracturing and weathering of the geological formations. The data are too few to enable a conventional modeling approach of this volcanic system. The implementation of the GLUE method in a numerical groundwater flow model allowed developing a stochastic analysis of the spatial distribution of the hydraulic conductivity and the recharge modalities of this complex volcanic system. The hydraulic conductivities range from 10⁻⁶ to 10⁻⁸ m∙s⁻¹ for the basalt and the rhyolite formations (values are yet generally lower for rhyolites) and are higher than 5 × 10⁻⁷ for the sedimentary formations. In addition, considering diffuse recharge as the main mechanism by which the precipitation reaches the aquifer results in more consistent groundwater head simulations than considering only indirect recharge. The average recharge amount estimated for the Goda aquifer system is 28 mm∙yr⁻¹. The results led to a numerical representation of this system, with the least uncertainty. This model was able to estimate the available water resources of this system. This result is important because the Goda system supplies water to the city of Tadjourah. Assessment of available resources is vital for the future development of this city. From a methodological point of view, the GLUE method proved very promising for water resources assessment in complex hydrogeological systems for which little data are available.

Keywords:

Water Resources, Stochastic Modeling, GLUE, Volcanic Systems, Goda Aquifer, Djibouti

1. Introduction

Understanding the hydrodynamic functioning of volcanic aquifer systems is confronted with the difficulty of dealing with their geometrical complexity and the high spatial heterogeneity of the recharge processes and flow parameters distributions [1] [2] [3] [4] . Volcanic aquifers have been studied quite extensively in developed countries [5] [6] , but few studies have been undertaken in developing countries [7] , where they are very often the only available water resources [8] [9] .

Given their strong heterogeneity, studying these hydrogeological systems require a high density and quality of data, which explains why such studies have been mainly carried out for systems equipped with well-developed monitoring networks. Volcanic aquifers for which access to data (in time and space) is limited are generally poorly described whereas they can be a crucial resource, especially in arid environments. The lack of data constitutes a real difficulty in assessing the resources of these aquifers [10] [11] .

More particularly, optimization procedures, which consist in minimizing the differences between simulated and observed loads, generally implemented in order to identify the parameters of the hydrodynamic models lead in some cases to the instability or non-uniqueness of the solution [12] [13] [14] [15] [16] . Strategies to better account for uncertainties should thus be implemented such as the GLUE (Generalized Likelihood Uncertainty Estimation) method [17] . This stochastic method consists in testing a large number of parameter sets and keeping all those which allow to simulate the observed hydraulic heads taking account of their uncertainties.

The GLUE procedure is quite recent and has been very little applied in Hydrogeology [18] [19] [20] [21] . However, this procedure would offer very interesting perspectives when aquifers are very complex with very little data available. This situation is frequently encountered in developing countries concerning volcanic aquifers. Under these circumstances, the little data available does not permit to conduct a classical and effective modeling of these systems. Yet modeling these systems is imperative, since they represent the only resources available, which must necessarily be exploited within a sustainable framework.

In this work, the GLUE approach was applied to the modeling of the Goda volcanic aquifer in the Republic of Djibouti (Horn of Africa).

2. The Study Area

2.1. Geomorphological Features

The Goda massif (Figure 1) is located to the west of the city of Tadjourah. This steep area consists of hills and narrow valleys. It is bordered to the south-east by the littoral and is limited to the north by the plateau of Dalha at the village of Adaylou. On the northern and western slopes of Day Mountain, which peaks at 1800 m, the altitude of the massif decreases steadily. The massif is limited to the south and west by Lake Assal (−150 m asl) and the fault zone of Makarassou. The central slope of the massif is home to Day’s primary forest and consists of three exoreic watersheds (Darriyou, Magalé, Asleyi). The relief is particularly steep with dips reaching more than 35˚ on the descent of Day Mountain.

Precipitations over the massif of Goda are higher than the average over the territory of Djibouti (150 mm/year). Indeed, more than 350 mm/year of rainfall are recorded on the high zones. The precipitation decreases until reaching 200 mm/year on the coast near the city of Tadjourah (Figure 1).

2.2. Geology

The Goda massif is located on the terminal part of the Danakil block bounded to the north in Eritrea by Precambrian sedimentary rocks [22] [23] . The tectono-volcanic activities following the opening of the Red Sea [24] [25] [26] set up

the volcanic rocks of the Goda massif. The oldest geological unit is the formation of the Mabla rhyolites dating of the Miocene [27] . This unit outcrops on the central and northern part of the massif. Normal faults of N160 or N50 direction intersect the series of Mabla rhyolites.

The Mabla rhyolites are covered to the west by Dalha basalts (4 - 9 Ma) and stratoid basalts (1.25 - 3.5 My) [28] . Dalha basalts are weakly faulted and subhorizontal whereas the stratoid basalt series is affected by normal N-S main faults. There are some outcrops of the Gulf basalts (1.25 - 3 My), of the Asal series basalts [29] and of the Mabla basalts, which are a contemporary series of the Mabla rhyolites. Sedimentary formations (<2 My), composed of detrital materials, are mainly located on the coast, but can also be intercalated in volcanic formations (Figure 2).

The aquifer system of the Goda massif is defined by three hydrogeological units consisting of Dalha basalts, Mabla rhyolites and sedimentary formations. On the top of the massif, in the village of Day, two boreholes drilled until depths of 323 and 358 m allowed to estimate their thickness (locally) to 300 m. The total thickness of the Mablas rhyolite is not known. Electric soundings performed at the mouths of the Darriyou and Magalé rivers have shown that the total thickness of the sediments could exceed 300 m.

2.3. Data

Over the whole Goda massif, 38 water points are available with piezometric data. They include 19 wells, 10 boreholes and 9 springs (Table 1). Pumping tests were carried out on only 5 boreholes, 4 of which are located in the Dalha basalts and only one in the sedimentary formations (Table 2). The water points are not evenly distributed over the study area. They are mainly concentrated in the central part of the massif and along the coast.

3. Methodological Approach

In order to evaluate the dominant recharge processes in the south-eastern slope of the Goda massif, a numerical modeling of the groundwater flow was implemented by considering two distinct scenarios: i) a diffuse recharge over the entire zone; and ii) an indirect recharge from the rivers beds. The MODFLOW- USG code [30] was used to solve the equations governing underground flows.

3.1. Geometry of the Volcanic System and Boundary Conditions

The Goda aquifer system is conceptualized with two layers: i) a top layer corresponding to a relatively productive zone with a thickness of 300 m; and ii) a semi-permeable lower layer that corresponds to a relatively unproductive zone in relation with the progressive clogging of deep fractures with a thickness of 1000 m. The top of the upper layer is defined by the DEM (Figure 1).

The model domain is limited to the south-east by the coastline where a fixed head of 0 m is set, to the south-west by Lake Assal where a fixed head of -150 m is set; to the northwest by a head contour line of 800 m, 30 km northwest from the summit of the Goda, where a fixed head of 800 m is set.

The area of the domain is 1800 km² and a regular discretization with a mesh size of 500 m × 500 m (11,600 mesh) has been considered. The boundary conditions are shown in Figure 3.

3.2. Hydraulics Parameters

A constant hydraulic conductivity (K) of 10⁻⁸ m∙s⁻¹ is assigned to the lower layer, considered to be semi-permeable. Lower values of K do not significantly modify the simulated heads, but increase the calculation times. For the upper layer, three zones of hydraulic conductivity are defined on the basis of the geological map (basalts, rhyolites and sedimentary formations) (Figure 4).

In the first modeling scenario, a diffuse recharge is considered according to three zones defined on the basis of the isohyets (>350 mm∙y⁻¹; 250 to 350 mm∙y⁻¹; >250 mm∙y⁻¹) (Figure 5). In the second scenario, recharge is applied only to the model cells corresponding to the drainage network. The three zones defined based on the isohyets are associated with different degrees of the network development and with distinct indirect recharge values (Figure 6). To prevent the accumulation of water above the topographic surface, mesh drainage at the level of the DEM is carried out where necessary.

3.3. Identification of Parameters

To identify the values of both hydraulic conductivity and recharge on the different zones defined above, a GLUE (Generalized Likelihood Uncertainty Estimation) approach is implemented. This method considers that several sets of input parameters can result in a concordance between the simulated and observed values of a given variable. The strategies adopted are the following:

1) 10,000 sets of parameters (hydraulic conductivities ranged from 10⁻⁴ to 10⁻⁸ m∙s⁻¹ and recharges ranged from 0 to 400 mm∙y⁻¹ for each of the different zones

defined in Section (3.2) are generated randomly according to a Uniform distribution;

2) a model performance indicator is defined, based on the Root Mean Square Error (RMSE) between the measured and simulated hydraulic heads on the 39 control wells;

3) 10,000 simulations (with each of the parameter sets) are run and the 250 sets of parameters resulting in the lowest RMSEs are kept to establish probability densities.

4. Results

4.1. Scenario 1: Modeling with Diffuse Recharge

The differences between observed and simulated heads for the 250 selected simulations vary significantly according to the control wells. These differences must be compared with the uncertainty in the topography, in particular due to the change in DEM resolution from 30 m to 500 m in modeling. Indeed, in zones with high vertical differences, this uncertainty reaches ± 100 m (Figure 7). The probability densities obtained for K indicate that volcanic formations are

less permeable than sedimentary formations. The K values obtained for the Dalha basalts range from 10⁻⁸ m∙s⁻¹ to 10⁻⁶ m∙s⁻¹ with a modal value of 10⁻⁷ m∙s⁻¹. The K values obtained for Mabla rhyolites are generally lower, with a modal value of 10⁻⁸ m∙s⁻¹. The K values obtained for sedimentary formations are higher than 5 × 10⁻⁷ m∙s⁻¹ and two modal values are distinguished, 10⁻⁶ and 5 × 10⁻⁴ m∙s⁻¹ (Figure 8). The probability densities obtained for the diffuse recharge indicate that this parameter is lower for the upper zone (values ranged from 10 to 30 mm∙y⁻¹; modal value of 20 mm∙y⁻¹) than for the middle zone (values higher than 20 mm∙y⁻¹; modal value of 40 mm∙y⁻¹). For the low zone, not any clear trend can be identified (Figure 9).

4.2. Scenario 2: Modeling with Indirect Recharge

Overall, the same trend as in the first scenario is observed for the second scenario, but on average the hydraulic conductivity values are slightly lower than the previous ones. Indeed, values between 3.4 and 8.5 × 10⁻⁸ m∙s⁻¹ for volcanic formations and a value of 3.4 × 10⁻⁶ m∙s⁻¹ for sedimentary formations are obtained (Figure 10). For this configuration, since the recharge is assigned only to the drainage network, the values required to adjust simulation and observation are

high. The upper zone receives 52.4 mm∙y⁻¹ and the middle and low zones 123.7 and 198.2 mm∙mm∙y⁻¹ respectively (Figure 11).

5. Discussion

The probability densities of the hydraulic conductivities of volcanic formations have well-constrained unimodal distributions. The values of hydraulic conductivity obtained for these formations remain on average lower than the values found for the rhyolites. However, for the indirect recharge scenario, the model displays for the basalts a hydraulic conductivity higher than one order of magnitude compared to the modal value obtained for the Mabla rhyolites.

The probability density distributions for sedimentary formations are less constrained and show two subsets of hydraulic conductivity fields for the two recharge configurations. Indeed, according to Jalludin and Razack [3] , the hydraulic conductivity of the sedimentary formations presents a great variation due to the difference in the constitutive materials of the alluvial deposits.

The only study to estimate recharge, based on field experiments (differential gauges), dates from 1982 and concerns the coastal aquifer of the Gulf basalts [31] . According to this study, this aquifer would receive only 5% of the annual raw precipitation. The amount of recharge estimated for the different zones depends entirely on the mechanism governing the recharge. Indirect recharge is the dominant mechanism in arid contexts [32] as for Djibouti volcanic aquifers [31] [8] , but chemical and isotopic analyzes have shown that aquifers could also benefit of diffuse recharge [33] .

Analysis of the different scenarios which have been tested, pointed out that the simulated hydraulic heads are closer to the observed heads in the case of a diffuse recharge scenario (Figure 12). The recharge of the Goda volcanic aquifer, corresponding to this optimal scenario, amounts to 50.8 × 10⁶ m³/year representing an average 28 mm/year over the whole area. This estimation of the recharge is important for the concerned area as it allows a better understanding of the resources that are available and their rational exploitation in the future.

6. Conclusions

The classical modeling approaches to estimate the hydrodynamic parameters and the recharge of strongly heterogeneous aquifers rove insufficient in a context of lack of knowledge and lack of representative data of the studied aquifer system. This is all the more true as the complexity of the concerned aquifer is high.

By testing a wide range of input data (hydraulic conductivity, recharge) for the different main configurations by which the water infiltrates towards deep aquifers, the GLUE approach made it possible to understand that volcanic formations (basalts, rhyolites) have rather low hydraulic conductivity, which order of magnitude is similar for all formations. By cons, the sedimentary formations present significant heterogeneities in hydraulic conductivity probably due to their composition.

The estimated recharge rate varies between 6% and 12% of the considered isohyet when diffuse recharge is the dominant mechanism. On the whole, the simulated heads are as close as possible to the observed heads when the precipitation infiltrates diffusely into the aquifer.

This modeling approach has proved to be relevant for a better understanding of the functioning of the complex volcanic system of the Goda rhyolitic massif. An assessment of water resources could be deduced. This result demonstrates the actual performances of the GLUE method in Hydrogeology. This work opens up promising prospects for the modeling of very complex aquifer systems, for which there are very little available data.

Cite this paper

Ahmed, I.M., Le Coz, M., Jalludin, M., Sardini, P. and Razack, M. (2018) Assessment of Groundwater Resources in a Complex Volcanic Reservoir with Limited Data Sets in a Semi-Arid Context Using a Novel Stochastic Approach. The Goda Volcanic Massif, Republic of Djibouti. Journal of Water Resource and Protection, 10, 106-120. https://doi.org/10.4236/jwarp.2018.101007

References