ABSTRACT

The present work applied the HydroGeoSphere (HGS) model in humid tropical area to Koué watershed scale to simulate flows in porous and fractured area of crystalline aquifers. It integrates rainfall, physiographic, fractures, hydraulic drilling and hydrodynamic data. The simulation of flows in porous area concerned 5 test zones. The input database of the model is implemented on a triangular grid in porous area using Gridbuilder software and interactive block grid in fractured area. In order to use the model in these two environments, boundary condition was set. The infiltrations rate of the earth layers is estimated in the order to 10⁻⁵ ms⁻¹. The model simulates the pumping with a good reproductivity of the drawdown profiles of groundwater at the drillings. The storage coefficients vary between 9.9 × 10⁻⁴ and 2 × 10⁻³. The hydraulic conductivities vary from 8.5 × 10⁻⁶ to 2 × 10⁻⁵. 73.9% of the drillings studied has a high hydraulic conductivity and shows a strong drawdown of the groundwater table. The study of the static levels of the ground water allowed indicating the distribution of the water resources in the drillings: 57% are deep in the first 10 meters, 36% between 10 and 20 m, and 7% in the higher level to 20 m deep in the earth.

Keywords:

Crystalline Aquifer, Hydrogeological Modelling, HydroGeoSphere, Water Resource, Koué Watershed, Côte d’Ivoire

1. Introduction

The water sector includes direct use and consumption of the resource, land drainage and irrigation, flood assistance, fisheries, use for industry, protection of environment, collection and treatment of wastewater and industrial effluents etc. The groundwater is one of the major themes of reflection on sustainable water resource management. In recent decades, many countries have faced many challenges in the availability and exploitation of water resources (high population growth, rapid urbanization, economic pressures and climate change). In such a situation, groundwater resources with low vulnerability to drought seem to be the best alternative to respond to this emergency and adapt to impacts (Koïta et al., 2018) . Assessment and integrated management of water resources is part of the National Development Plan (NDP) and this management is done by hydrologic and hydrogeologic watershed.

In Côte d’Ivoire, apart from the coastal sedimentary basin, the main groundwater resources are found in the fissured rocks of the Precambrian basement and in the alterites already exploited by most village wells (Biémi, 1992) . The semi-mountainous region studied, consisting mainly of crystalline rocks, has been the subject of several studies in geology (Azuelos et al., 1978; Camil, 1984; Kouamélan, 1996) .

Thanks to the projects “Fouta Djalon” and “Radar-Water”, the cartography of the geological structures affecting the Archean rocks was carried out as part of the evaluation and management of water resources (Kouamé, 1999; Lasm, 2000; Assoma et al., 2013) . These studies have located the high hydraulic potential areas for precise hydraulic drilling implementations at high flow.

Estimation of hydrodynamic parameters has been followed by methods of application, including logging and geoelectric measurements (Kelly, 1977; Perdomo et al., 2014) . It is the same for the combined estimate of the specific performance and natural recharge in an unconfined aquifer (Maréchal et al., 2006) . Developments in Magnetic Resonance Sounding (MRS) methods particularly allow a description of aquifers and provide direct links with hydrogeological properties (Legchenko et al., 2002; Vouillamoz et al., 2014) . Authors (Koïta et al., 2018) have combined these magnetic resonance methods with ground level fluctuation data for estimating groundwater storage change.

The recent study (Dewandel et al., 2012; Rezaei & Mohammadi, 2017) proposed methodologies for regionalizing specific yield at the watershed scale which combine at cell scale (pumping scale) water-table fluctuation and groundwater-budget techniques in the absence of recharge from rainfall. The particular and innovative character of the HydroGeoSphere model applied to this study is that it integrates both rainfall, physiographic and hydrodynamic data for the estimation of groundwater resources, the characterization and hydrogeological modeling of aquifers.

In fractured area, fluid circulation and material transport are strongly dependent on the geometry of the fracture network. To contribute to advanced knowledge of aquifer systems, the use of mathematical models has become unavoidable. These models incorporate multi-source data to optimize water resources management. These digital models provide relevant tools for analysis and simulation. This solves the problems of underground flow and contaminant transport in highly heterogeneous and complex geometry systems (2D and 3D). The quantification of groundwater resources and groundwater recharge (infiltration or exchange between groundwater and surface water) inside the watershed is necessary for efficient management of this resource.

This study aims to ensure a better understanding of hydrogeological functioning at the Koué watershed using the HydroGeoSphere model, in order to contribute to the rational management of groundwater resources. It proposes simulations on five test sites at the watershed scale of 2045 km². The methodological approach is based on a spatial discretization of the input data of the HGS model and the implementation of which requires the definition of the boundary conditions.

2. Study Area

The Koué River is one of the tributaries of the Sassandra River. The Koué watershed (2045 km²) studied is located in the semi-mountainous region of Man in western Côte d’Ivoire (Figure 1) and has a very contrasting and very high relief, which rises to 1278 km². It consists essentially of a part of the massif of the chain of Touras mountain (500 to 1200 m), dotted with means and high plateaus and plain (Assoma et al., 2013) . Annual rains vary from 500 to 2200 mm with an average of 1467 mm of the period 1976 to 2002. The geological formations consist of Precambrian rocks and are mostly crystalline or crystallophyllous in nature composed of gneisses, migmatites, granulites, charnockites, magnetite quartzites and associated formations.

Figure 2(a) shows an example of a geological section made from village

hydraulic techniques and aquifer models in the basement of the study area. The aquifer model of the Koué watershed (Figure 2(b), Figure 2(c)) is a bilayer system formed by an aquifer of alterites (A) overcoming the fracture aquifer (B). Hydrodynamically, the capacitive and draining functions coexist within these two reservoirs. The altered area is essentially capacitive and the fractured area has a strong conductive tendency (Kouamé, 1999) .

3. Data and Materiel

Integrated water resources simulation required soil textural classes from the units of pedology (Perraud & Souchère, 1969) , digital terrain model (DTM) and the hydraulic drilling data. These Hydraulic data come from the datasheets of 103 drillings (1981 to 1982) and consists of thicknesses of alterites (0 to 55 m), thicknesses of the soil layers, operating flows (m³h⁻¹), static levels of groundwater (m) and total depths of structures (m).

The hydrodynamic parameters studied consist of the porosity (θs) which is the saturation water content, the hydraulic conductivity Kf (ms⁻¹), the fracture induced permeability (Kouamé, 1999) , the storage coefficient (S), the residual water saturation (θr) and the parameters α and β of Van Genuchten. These parameters come mainly from the study of Carsel & Parrish (1988) due the similarity with the soils of the Koué watershed.

The daily rainfall (1980 to 2000) at the Man station comes from the Development Company and Airport, Aeronautic and Meteorologic Operating (SODEXAM). The rainfall data were used to calculate effective rains and thus average daily infiltration and recharge of the water table and groundwater. The software used consists of:

- QGIS-2.18 for the mapping of soil types;

- ENVI-5.1 for DEM extraction and export in Gridbuilder accessible format;

- Gridbuilder has been used to mesh and convert maps (DEM, flows of drilling exploitation, thickness of alterites, total depth of structures) in a accessible format to HGS and Tecplot models;

- HGS model (grok, hsplot, hydrogeosphere) for the simulation of surface and underground flows;

- The Tecplot model for viewing and displaying the simulation results.

The HGS model provides a capability of rigorous simulation that combines hydrology—water quality—subsurface flow and transport. It is fully integrated with a set of user interface tools. It is a simulator (in 3D) of finite elements of control volume designed to simulate the entire terrestrial part of the hydrologic cycle.

4. Methods

4.1. HydroGeosphere Model Input Database

4.1.1. Discretization of Physiographic Data

HGS is a finite element model based on the discretization of the space of the domain where the flow is studied. A model of finite element solves a matrix system based on flow equations for all of grid nodes. Two types of grid have been used for the simulation of the different lithological layers (Figure 3):

- a triangular grid in 2D (x, y) for the porous area performed using the GridBuilder Tool (McLaren, 2008) . The grid is made on the basis of the textural classes of the soil. Gridbuilder does not allow to choose the desired size of elements but generates a default grid that can be refined;

- an interactive block of grid in 3D (x, y, z) in a fractured environment, suitable for breaks or wells in which fine grid are required (Therrien et al., 2010) . This type of grid can be defined by the user as desired in the three main directions.

4.1.2. Discretization of the Soil Texture Profile

The textural classes of the soil for the superficial layers (Figure 4) are edited from the soil map. Five soil units have been identified (Table 1) and take into account the particle size of soils that influences hydrologic processes (infiltration and surface flow).

For the use of these soils in the HGS model, an analysis according to the average texture composition (sand, silt and clay) was carried out. It refers to the textural triangle nomenclature USDA (United States Department of Agriculture). This approach identified four (4) textural classes of soil that correspond to the upper layers of the basin: silty sands (15.1%), sandy loams (19.7%), loamy-clay

sands (40.7%) and clay sands (24.5%). The simulations are performed separately in each of these four textural classes.

The identification of the profile of the four textural classes is based on the data of the hydraulic drilling datasheets. For each textural class within the Koué watershed, several descriptions of the drilling profile were used to situate the different thicknesses and the names of the identified lithology (Table 2):

- three lithologies in the area of Facobly, Flansobly and Siambly;

- two lithologies in the area of Mangouen and Gbonné;

- clay is almost the second layer of the different textural classes of the soil, except in Gboné where she is absent.

4.1.3. Discretization of Hydrodynamic Data

The hydrodynamic parameters used in this study are also used by many models such are KINEROS (Woolhiser et al., 1990) , STICS (Brisson, 1989) , HYDROTEL (Fortin et al., 2001) and HydroGeoSphere (Therrien et al., 2010) to estimate the apparent physical parameters and water parameters of soils. The model is numerically constructed according to the following parameters which are declined in spatial or algebraic discretization:

- parameter of “soil”: the correspondence between textural class and hydrodynamic properties is related to the study of Rawls & Brakensiek (1989) and Carsel & Parrish (1988) . According to Rawls et al. (1982) , these results can be used to study theoretic models by comparison with a large set of experimental data;

- parameter of “fractured environment”: the drilling datasheets consulted indicate pumping of short times (8 to 12 hours) but no long times (24 to 72 hours). The calculated hydrodynamic parameters relate more to the structure and the action of the solicited fractures during the pumping (Kouamé, 1999) . As the watershed is intensely fractured, the underlying saprolite (soft part of the weathering profile) can be assimilated to an equivalent and continuous porous environment for the calculation of hydrodynamic parameters;

- parameter of “storage coefficient (S)”: this parameter could not be evaluated for lack of observation piezometers in the study area. In the absence of any rigorous determination from isolated drilling, the coefficient S values between 10⁻⁴ and 10⁻³ are suggested for the reservoirs of fractures (Faillat, 1986) ;

- parameter of “hydraulic conductivity (K_f)”: in a space system defined by a disc of diameter D, containing a set of mega-fractures (Figure 5), the hydraulic conductivity (K_f) of each fracture is calculated on the basis of Franciss method (Equation (1)).

$K_f\left({\text{ms}}^{-1}\right)=\frac{e}{D}$ (1)

with, K_f: hydraulic conductivity of the mashed area (ms⁻¹); e: thickness of the mashed area (m); D: diameter of the disc (m).

Assuming the fracture is continuous inside the disc and the thickness (e) of the mashed zone is proportional to the fracture length (L), the hydraulic conductivity (K_f) is derived from Equation (2).

$K_f=\frac{CL}{D}$ (2)

with, L(m): length of the fracture; C = e/L: empirical proportionality coefficient between the thickness e(m) of the mashed area and the length L of a mega-fracture.

In Côte d’Ivoire, this coefficient of proportionality (C) varies from 4.3 × 10⁻⁴ to 1.2 × 10⁻² with an average of 2.8 × 10⁻³ (Kouamé, 1999) . Figure 6 presents the mean induced permeability by fractures and the values of the hydraulic conductivity.

The conductive nature of the porous area is appreciated by the experimental law of Darcy which governs the flow of an incompressible fluid in porous earth. It leads to define the permeability K of the considered ground in relation to the phenomenon of infiltration.

The infiltration rates obtained are calculated using HGS model according to Equation (3).

$V=\frac{P_L}{t}$ (3)

with, P_L(m): piezometric or dynamic level of groundwater table; t: simulation times (s).

4.2. Implementation of the HydroGeoSphere Model

4.2.1. Boundary Conditions and Simulation in Porous Area

For the pedological units taken into account, the contours of the domain were chosen so as to integrate the studied boundary conditions:

- the first boundary condition sets a hydraulic charge (0 m) with the assumption that the underground flow is zero (condition of Neumann);

- the second boundary condition is the imposed flow: This is effective rainfall (0.3 mm per day) representing the mean of daily infiltration intended for the recharge of the groundwater and deep aquifers. It is estimated from the water balance with an average annual infiltration of 110.2 mm.

To model the groundwater flow in saturated conditions, HGS model considers the following assumptions:

- the fluid is essentially incompressible;

- the porous area and fractures (or macro-pores), if they exist, are non-deformable;

- the system is according to isothermal conditions;

- the air phase is infinitely mobile.

Table 3 presents the hydrodynamic parameters of the different layers (on the basis of the textural class and the datasheets of the drillings).

For the study of flows and infiltrations through the hydrogeological system, 80 hydraulic drilling were used. Some of these parameters (porosity, saturation, α and β) are taken from the study of Carsel & Parrish (1988) .

The use of the HydroGeoSphere model in a porous area requires the water retention profiles of the aquifer (degree of saturation in terms of pressure) and the relative hydraulic conductivity versus saturation of the different materials (Figure 7). This integrates the parameters α and β of Van Genuchten. The simulation is performed in the beginning of the rainy season (march).

4.2.2. Boundary Conditions and Simulation in a Fractured Area

The simulation in fractured area (Figure 8), considered as an “Equivalent Porous Area” (Nastev et al., 2004) is type of axisymmetric. It is made out within a

radius of 1000 meters around the drilling. The use of the datasheets of the hydraulic drilling made it possible to have the following parameters:

- measures of drawdown of the groundwater table. They are provided by pumping tests, which are the short well tests (6 h), carried out on 23 isolated drillings in the areas of Facobly and Gbonné;

- alteration thicknesses, bedrock depth and air-lift flow.

One of the boundary conditions consist in determining the static levels of the groundwater from the drilling datasheets due to the lack of piezometers. The pumping is done from a vertical well, located at the coordinates (X = 0, Y = 0, Z = drilling depth). Each drilling has a radius of 0.082 m with a characteristic pumping flow. Along the OX axis, the triangulation elements are refined using the HGS model in the area near the well. The initial size of the elements is 0.01 m and then increases to a maximum of 100 m. This approach is based on the method of Theis (1935) in transient flow as Equation (4), for calculating the drawdown (S_c) in time of descent of a pumping of short duration.

$S_c=\frac{Q}{4\pi T}W\left(u\right)$ (4)

with, $W\left(u\right)={{\displaystyle \int }}^{\text{}}\frac{{\text{e}}^{-v}\cdot \text{d}v}{v}$ and $u=\frac{r^2\cdot S}{4Tt}$ ; Sc: folding in the piezometer (semi-captive

sheet); Q: pumping rate of the well; T: transmissivity; S: storage coefficient; t: time (s); W(u): can be calculated from a table.

It expresses the drawdown of groundwater table at a distance (r) from the pumping well with a constant flow rate Q after a pumping time (t) and thus respects the followings conditions of applicability:

- the aquifer must have an infinite lateral extension;

- the diameter of the well is negligible;

- the well is perfect;

- the flow to the well is transient;

- the aquifer is active, homogeneous, isotropic, with a constant thickness over the entire area influenced by the pumping.

Fracturing at all scales becomes the dominant factor that determines the amount of water that can be contained in the rock and possible flows. In practice, the storage coefficient estimated from the sampling flow per drilling is used. The procedure consisted in comparing the theoretical drawdown profile performed using HGS model and the experimental profile versus time, using the values of the storage coefficient (S) and permeability (K), because HGS model does not take into account the transmissivity (T). The selected values of (S) must vary between 10⁻⁴ and 10⁻³ in order to be consistent with the values of fissured environments in the region (Lasm, 2000) . The selected values of hydraulic conductivity (K_f) are based on the study of Kouamé (1999) .

The choice of the two test area (Facobly and Gbonné) for the simulation is justified because these zones are favourable to strong potentiality of water affected by bigger fairways couloirs of strong intensity of flow (Assoma et al., 2013) . This is conducive to the process of storing a significant amount of water resource in highly mashed or altered areas. The channelization of water sometimes passes through several high-flow areas. This shows that the aquifers are continuously supplied and interconnected by through major fractures.

5. Results

5.1. Simulation of Flows in Saturated Porous Area

The simulation in permanent scheme of flow in the saturated zone allowed interpreting the temporal evolution of the dynamic level of the groundwater in each of the test zones of the watershed. Figure 9 shows the case of modelling the saturated flow zone in a three-layer of hydrogeological system of Facobly test area.

The temporal distribution of the dynamic level of the groundwater, assimilated to piezometric level during readjustment, is unevenly distributed in the lithological layers. These dynamic levels simulated of the groundwater tables allowed to characterize the infiltration rate of water (Table 4), which represent an indicative values of the permeability of each lithological layers.

The average infiltration rate obtained is 4.6 × 10⁻⁵ ms⁻¹. In porous area, infiltration rates are twice as high in the Gbonné area. In the rest of the zones (Siambly, Facobly, Flansobly and Mangouen) the low infiltration rates reflect slow groundwater recharge. The correlation between the lithology and its infiltration rate could be established (Table 5).

The relative hydraulic conductivity (degree of certainty) depends on the degree of saturation (following the rainwater infiltration) which is closely related to the pressure of the hydraulic charge of the groundwater (Figure 10):

- hydraulic conductivity is higher in sand, followed by clay sand and very low in clays because they retain more water;

- the pressure of hydraulic charge is higher in clay, lower in clay sand and the sand because they retain less water;

- considering the infiltration rate values of 4 × 10⁻⁵ ms⁻¹, the water supply to the aquifer of the Koué watershed catchment area is declared good;

- the values of hydraulic conductivity range from 8.5 × 10⁻⁶ ms⁻¹ to 2 × 10⁻⁵ ms⁻¹ with an average of 1.4 × 10⁻⁵ ms⁻¹. This expresses a fairly good performance of the HGS model and its applicability in humid tropical environments.

5.2. Simulation of Flows in Fractured Environments

The simulation is performed in permanent scheme in the saturated zone. This is an axisymmetric simulation (2D) that generated 3534 nodes in total and 1680 elements (grid) for each case in fractured environment. The simulations performed at the first levels of the drawdowns are quite satisfactory for all of 23 drillings (Figure 11).

The storage coefficients (S) of the fissured aquifers studied vary between 9.9 × 10⁻⁴ and 2 × 10⁻³ with an average of 1.5 × 10⁻³. The distribution of the hydraulic conductivities of the 23 hydraulic drillings according to their intensity provided four classes values (Table 6). In this distribution, 73.9% of the drillings studied offer high hydraulic conductivity and express a strong drawdown.

5.3. Water Resource of Aquifers in the Bedrock Environment

The map of static levels of the aquifer shows the distribution of groundwater resources, based on a study of 103 hydraulic drillings in the Koué watershed (Figure 12). These groundwater resources are distributed according to depth: we take note that 57% are contained within the first 10 meters, 36% between 10 and 20 meters, 7% beyond 20 meters.

The presence of an altered unsaturated cover is an unfavourable factor for the recharge of the fractured aquifer. Indeed, if the hydrostatic level is located in the fractured area, the altered unsaturated layer can immobilize the fraction of rain effective to fill its moisture deficit. The renewal of water resources in this case is not ensured.

6. Discussion

Knowledge of the characteristics of aquifers and the conditions of their recharge are useful for the assessment and management of water resources in order to meet the water needs of populations. In the context of groundwater management or for hydrogeological modeling, precise estimates of the spatial variation of hydrogeological properties, in particular effective hydraulic porosity and conductivity

as advocated by Dewandel et al. (2017) , on two large crystalline aquifers in the watersheds of southern India. These studies are based on the method of regionalization of 3D effective porosity.

In this study, simulation within the porous area allowed to determine infiltration rates of the soil layers (silty sand, sandy silts, sandy clay and clay). In the soils of the Koué watershed studied, the infiltration rates of sand, silty sand and sandy silt estimated correspond to the proposed standard (Table 7). However,

there are small differences in the observed values due to the experimental protocol used to set the values of the standard tests.

In addition, the soils studied are deep and the infiltration rate was determined from the effective rainfall. Indeed, the tests must be carried out on shallow soils between 30 and 80 cm deep, in normal meteorological conditions, without rain or frost with clear water, free of a high concentration of organic products and preferably at night. The hydraulic conductivity values resulting from the simulation are in the same order of magnitude as those estimated by Oularé et al. (2016) in the N’Zo watershed in the Man region with estimated values ranging from 1.1 × 10⁻⁶ and 2.4 × 10⁻⁵ ms⁻¹. These authors have adopted a reverse modelling approach, the principle of which consists in calculating the parameters of the aquifer that are difficult to access in the field based on observations that can be easily collected in the field.

Simulations of pumping in fractured environments using HGS model show a good estimate of drawdowns compared to those measured in the drillings of the Facobly and Gbonné areas. The storage coefficient varies between 9.9 × 10⁻⁴ and 2 × 10⁻³. The hydraulic conductivity is between 8.5 × 10⁻⁶ and 2 × 10⁻⁵. These results are in agreement with those of de Faillat (1986) for whom the drawdown values are between 9 × 10⁻⁵ and 1.1 × 10⁻³. It is the same for the hydraulic conductivity values of Kouamé (1999) who established an average induced permeability map of the western mountainous region of Côte d’Ivoire.

Based on the infiltration rates calculated using the HGS model, the conductive nature of the porous area and the equivalent porous area (fractured) through an imposed flow could be demonstrated in this study. The process of infiltration of the fraction of rain for continuous recharge of groundwater flow and thus groundwater recharge is a determining factor in the assessment of groundwater resources.

Linear pressure drops (or hydraulic gradient) and the fractures opening contribute to a better knowing of the aquifer environment: this hydraulic gradient through the system (which is inversely proportional to the importance of capillary processes), is the parameter with the greatest impact on the flow leaving (Pouget, 1998) . In permanent of regime, fracture opening is the parameter that most influences flow, regardless of network fracture and saturation conditions (Pouget, 1998) .

The HydroGeoSphere model offers the advantage of using uniform, variable and interactive types of grids (Therrien et al., 2010) . Generally, in simulating unsaturated flow conditions, it is desirable to consider a finer grid size in order to better approach the solution and take into account variations in hydraulic conductivity and water content versus pressure (Farouk, 2009) .

7. Conclusion

The issue raised at this work was to contribute to a better knowledge of the hydrogeological functioning of the Koué watershed using the HydroGeoSphere code, which uses the finite element method. Thus, the integration of physiographic, hydrodynamic and rainfall data into the HGS model made it possible to simulate flows in porous and fractured environments in the Koué watershed. This model easily integrates all these data organized in a Relational Database Management System (RDBMS) and represented as tables in a GIS (Geographic Information System). The method adopted consisted of data discrediting and setting the input data from the HGS model. The implementation of this model is based on the definition of boundary conditions for simulations in porous and fissured environments.

The average water infiltration rates of the lithological layers calculated using the HGS model are 4.6 × 10⁻⁵ ms⁻¹ and close to those of the standard tests. Infiltration is more noticeable in the Gbonné area than in the rest of the Koué watershed. In general, the model simulates pumping well with good reproducibility of the drawdown profiles of each drilling. The average values of storage coefficient are 1.5 × 10⁻³ and that of the hydraulic conductivity is 1.43 × 10⁻⁵ ms⁻¹. Overall, 73.9% of the drillings studied offer a high hydraulic conductivity and reflect a high drawdown. The distribution of statistical levels of the groundwater table made it possible to indicate the distribution of water resources in the drillings versus depth.

Drinking water supply in the test areas studied is provided by groundwater resources contained mainly in the fractured basement. These results, while important, reflect the complexity of the fissured environment in water resource management at the watershed scale. A better understanding of hydrogeological functioning is fundamental to the management and protection of groundwater resources. Hydrogeological characterization is required to define the subsoil structure and identify the main processes that contribute to aquifer productivity and water quality. To this end, field observations and the use of geophysical, hydraulic and hydrochemical methods should be considered.

Acknowledgements

We sincerely thank the National Office of Drinking Water and the Direction of Human Hydraulics of Abidjan for the availability of data.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Assoma, V. T., Konan, M. Z., Adon, C. G., & Kouamé, F. K. (2019). Integrated Numerical Simulation of Water Resources in Crystalline Aquifers at Koué Watershed Scale (Western Côte d’Ivoire). Journal of Geoscience and Environment Protection, 7, 202-221. https://doi.org/10.4236/gep.2019.77014

References