Abstract
We made use of groundwater flow and mass transport equations to investigate the crucial potential risk of water pollution from hydraulic fracturing especially in the case of the Karoo system in South Africa. This paper shows that the upward migration of fluids will depend on the apertures of the cement cracks and fractures in the rock formation. The greater the apertures, the quicker the movement of the fluid. We presented a novel sampling method, which is the combination of the Monte Carlo and the Latin hypercube sampling. The method was used for uncertainties analysis of the apertures in the groundwater and mass transport equations. The study reveals that, in the case of the Karoo, fracking will only be successful if and only if the upward methane and fracking fluid migration can be controlled, for example, by plugging the entire fracked reservoir with cement.
1. Introduction
In the recent decade, shale gas has become one of the mainly functional natural gases for industrial countries. For instance in USA, it was supposed in 2009 that natural gas demand was accepted to augment from 23 tcf per year to 30–34 tcf per year in 2025 [1–5]. However, usual gases were not able to keep happy such a need. Consequently to gratify this demand, the alternative gas sources such as shale gas were expected to be the most important paraphernalia of this construction. It is perhaps important to recall that shale are fissile rocks composed of layers of fine-grained sediments [5–10]. Several techniques have been put in place to extract this shale gas from the fissile rocks. By using the advanced techniques of horizontal drilling and hydraulic fracturing, it now seems to be economically feasible to extract natural gas from the Marcellus shale [11]. Even though these techniques are well recognized, they are not without potential risk. Hydraulic fracturing uses high-pressure solutions to create and prop open fractures in rock to improve the flow of oil, gas, or water [12]. More than 750 different chemicals, ranging from benign to toxic, have been used in hydraulic fracturing solutions [12]. Although these additives are less than 2% by volume of the total fracturing fluid, hydraulic fracturing is a water-intensive process and at least 50 m3 of chemicals would be used for a typical 10,000 m3 hydraulic fracturing project [11, 12]. The crucial unknown is the potential risk of water contamination from hydraulic fracturing especially in the case of the Karoo system in South Africa.
2. Statement of the Problem
Water is not only the most abundant substance on earth, but also the substance on which all forms of life depend. No wonder then that man has always been preoccupied with precious resource. The largest volume of water on the earth is in the oceans (which cover nearly 75% of its surface), but this water is not potable or suitable for domestic and industrial purposes. Man, therefore, always had to rely on other fresh water sources to satisfy his needs for potable water.
The science on contamination of drinking water from shale gas drilling, fracking, and production, is recent and incomplete. A peer-reviewed, archival journal study from Duke University [13] found apparent migration of substantial amounts of methane from gas wells to private water wells as far out as 1000 m in the Marcellus play in Pennsylvania [12]. This is illustrated in the Figure 1 below. It is therefore important to investigate the possible effect of these pollutions to the shallow and deeper aquifer after the closure of the fracking area. The questions that arise at this level are the following: what happen after closure of the fracking zone if there is a leakage of the capped borehole? At which extent the private borehole of the neighbourhood will be affected? These questions can find answer only by presenting a suitable model for this situation. Figure 1 shows approximately the possible situation that takes place during the fracking process, in particular, the potential gas migration paths along a well. However, Figure 2 shows the insufficient cement coverage leading to possible migrations of gas along paths well.
To achieve this, we consider the case where 1000 ha was fracked from 30 wells as indicated in Figure 3.
The following assumptions will be considered in this study.(1)Measure the potentiometric measure at the bottom of the well (this will give the total piesometric head of the organic shale (i.e., hydrostatic + hydrodynamic pressures) which will include the artesian pressures of the Karoo).(2)Because the organic Ecca shale is over-pressurized (from Soekor wells), gas and water will flow from the well.(3)Now (a) if pressures rebuild 100%, the organic shale is not enclosed by impermeable boundaries and the fresh water in the area will be polluted sooner or later from the deep organic water. (b) But if pressures do not recover, the well is situated in a closed system with impermeable boundaries and there will be no water pollution from depth.(4)The over-pressured organic shale at, for example, 3000 m.(5)The potentiometric head is higher than the water-table in the fresh water aquifer.(6)The organic shale is not bounded by impervious boundaries.(7)In all the boreholes, the cement annulus will crack and form preferential flow paths given enough time and are closed on top.(8)There are private boreholes in the neighbourhood of the area.(9)The aquifer is confined, such that we have a permeability barrier, such that, after the closure of the boreholes, as the fluid pressure will increase in a rock, the fluid pressure approaches the lithostatic pressure and the forces act at the sediment grain contact.
3. Mathematical Formulation
The mathematical equation describing the flow of water via an aquifer can be found in [16–19]. In this paper, we use the simple analytical solution describing the relationship between the apertures and the discharge rate. It is assumed that the upward flow of water along the faulty cement annuli can be approximated by the well-known cubic law (parallel plate model for fractures). We can represent a fracture as a planar void with two flat parallel surfaces. The hydraulic conductivity of this fracture is defined as follows: where is the fracture aperture, is the density of water, is acceleration due to gravity, and is the viscosity of water.
The mean groundwater velocity through the fracture can be calculated as the product of the fracture hydraulic conductivity and the hydraulic gradient: where is the hydraulic gradient.
The transmissivity of an individual fracture is then and the flux along the fracture is where is flow in m3/d per m width.
The validity of the cubic law for laminar flow of fluids through open fractures consisting of parallel planar plates has been established over a wide range of conditions with apertures ranging down to a minimum of 0.2 μm. Artificially induced tension fractures and the laboratory setup used radial as well as straight flow geometries. Apertures ranged from 250 down to 4 μm, which was the minimum size that could be attained under a normal stress of 20 MPa. The cubic law was found to be valid whether the fracture surfaces were held open or were being closed under stress, and the results are not dependent on rock type. Permeability was uniquely defined by fracture aperture and was independent of the stress history used in these investigations. The apertures in this study are considered uncertain because, it is very difficult even in the field or real world problem to measure accurately the apertures. The next section is therefore devoted to the discussion underpinning the evaluation of uncertainties in this model.
4. Parameter Uncertainties Analysis
Parameter uncertainty can be defined as uncertainty that arises in selecting values for parameters in the various models. There are many parameters in this assessment that are uncertain. First, there are insufficient data about the site climatic, geological, and hydrological conditions. As a result, such parameters as sorption coefficients, moisture content, river flow rate, river depth and width, hydraulic gradient in the aquifer, and erosion rate are taken from the general literature. Some parameters used need to be specified more accurately, for example, evaporation or distance between the disposal facility and the river and between the disposal facility and residences. On the other hand, the sensitivity analysis aims at quantifying the individual contribution from each parameter’s uncertainty to the uncertainty of outputs. Correlations between parameters may also be inferred from sensitivity analysis. It is a frequent routine and recommended to perform the uncertainty and sensitivity analysis in tandem [20–23]. In this section, we present a discussion underpinning the parameter uncertainties analysis of the solution of the transmissivity, discharge rate, and velocity as function of aperture.
4.1. Samples Generation
The LHS method [24] is a type of stratified MC sampling [25]. The sampling region is partitioned into a specific manner by dividing the range of each component of . We will only consider the case where the components of are independent or can be transformed into an independent base. Moreover, the sample generation for correlated components with Gaussian distribution can be easily achieved [26]. As originally described, in the following manner, LHS operates to generate a sample size from the variables . The range of each variable is partitioned into nonoverlapping intervals on the basis of equal probability size 1/. One value from each interval is selected at random with respect to the probability density in the interval. The values thus obtained for are paired in a random manner with the values of . These pairs are combined in a random manner with the values of to form triplets, and so on, until a set of -tuples is formed. This set of -tuples is the Latin hypercube sample. Thus, for given values of and , there exist possible interval combinations for a LHS. A 10-run LHS for 3 normalized variables (range ) with the uniform p.d.f. is listed below. In this case, the equal probability spaced values are .
4.1.1. Efficiency of LHSMC
Consider the case that denotes an -vectors random variable with p.d.f. for . Let denote an objective function given by . Consider now the following class of estimators: where is an arbitrary known function and . If , that is, if is a fixed point for , then represents an estimator of . If , one obtains the th sample moment. By choosing is a step function), one achieves the empirical distribution function of at the point . Now consider the following theorems.
Theorem 1. If ’s are generated by LHS method, then the statistic in (5) is an unbiased estimator of the mean of , that is,
In this paper, we present a modified Latin hypercube sample called the Monte Carlo hypercube sampling method (MCHSM), and the method is presented in the next subsection.
4.2. Proposed Methodology (Combination of Monte Carlo and Hypercube Sampling)
It was demonstrated that the hypercube sample method was more efficient and less time consuming than the Monte Carlo simulation. However, this Monte Carlo simulation still presents some worth. In this section, we propose a methodology that combines both the Monte Carlo simulation and the Latin hypercube sampling as follows. Assuming that the uncertain parameter is and ranges within , then the first step in this method consists of generating the sampling via the Monte Carlo sampling within . The next step is to reduce the number of sampling by calculating the mean, the variance, and the standard deviation of the generated sample. These statistical parameters can then further be used to construct a distribution function, for instance, the normal distribution. With constructed distribution in hand, one can further apply the hypercube sample method to generate the final samples.
5. Applications
Iman and Conover [27] applied the LHS approach to cumulative distribution function (c.d.f.) estimation of the three computer models: environmental radionuclide movement, multicomponent aerosol dynamics, and salt dissolution in bedded salt formations. They reported a good agreement of c.d.f. estimations. In this section, the application of Monte Carlo Latin hypercube sampling to groundwater pollution will be discussed. In agreement with the real world problem, we assume that unknown parameters in (2), (3), and (4) are boundaries, that is, . Then according to the Monte Carlo Latin hypercube technique, we first generate sample via the Monte Carlo sampling and this is presented below in the histogram. We generated the sampling of apertures via the Monte Carlo and we represent it in Figure 4 below, where -axis represents the possible values of apertures. In Figure 5, we present the cumulative distribution function of apertures and their respective probabilities. Finally, in Figure 6, we present the Normal distribution of the generated apertures via Monte Carlo simulation.
According to the (MCHSM), we next generate a final sample of apertures from the cumulative distribution function. In this sample, we have generated 26 apertures using the cumulative function of the apertures generated via the Monte Carlo simulation, to each aperture, and we have associated a probability and the graphical representation is given below. Figure 7 shows the selected apertures obtained from the Latin hypercube sampling and of course their corresponding probability.
Using (2), (3), and (4) and expressing the relationship between velocity, transmissivity, and discharge rate, the values of the selected apertures can be used now to evaluate their correspondent transmissivity, velocity seepage, and discharge rate. The relations have been depicted in Figures 8, 9, and 10.
5.1. Cumulative Functions
The distribution for the transmissivity, velocity seepage, and discharge rate is presented as a cumulative distribution function (CDF) or as a complementary cumulative distribution function (CCDF), which is simply one minus the CDF. Hence, in our case the cumulative distribution function can be approximated as follows: where And is the probability that a value larger than will occur. The distribution function approximated above provides the most complete representation of the uncertainty in the transmissivity, velocity, or discharge rate that is derived from the distributions. Figures 11 and 12 show the cumulative functions of the transmissivity, discharge rate, and velocity.
Figures 13 and 14 show the normal distribution of the selected discharge rate and transmissivity.
5.2. Variance of the Sampling and Repeatability Uncertainty
5.2.1. Variance of the Sample
The form of estimator of the variance of is given by The goodness of an unbiased estimator can be measured by its variance. The variance approximated here provides a summary of this distribution but with the inevitable loss of resolution that occurs when the information is contained in 20 numbers [28].
5.2.2. Repeatability Uncertainty
It is important noting that repeatability uncertainty is equal to the standard deviation of the sample data [29]. In the case under investigation, the mathematical expression is given as follows:
5.3. Develop the Error Model
An error model is an algebraic expression that defines the total error in the value of a quantity in terms of all relevant measurement process or component errors. The error model for the quantity , that can be transmissivity or discharge rate, can be calculated with the formula below: where is the error in the transmissivity or discharge rate; is the error in measurement of viscosity of water; is the error in the measurement of aperture; is the error in measurement of hydraulic conductivity; is the error in the measurement of gradient; is the error in the measurement of density of water; and , and are the first-order sensitivity coefficients that determine the relative contribution of the errors in , and to the total error in . For this purpose, we chose the following definition of error: Then, Here, we also chose
5.4. Uncertainty in Quantities or Variables
The uncertainty in a quantity or variable is the square root of the variable’s mean square error or variance. In mathematical terms, this is expressed as follows: Providing that the correlation coefficients for the error in , and are equal to zero.
5.5. Skewness and Kurtosis Tests
Descriptive statistics, such as skewness and kurtosis, can provide relevant information about the normality of the data sample. Skewness is a measure of how symmetric the data distribution is about its mean. Kurtosis is a measure of the “peakedness” of the distribution [27–29]. In mathematical terms for the case under investigation, these are expressed as follows: since are our sampled functions from a sample of size 26, with mean and standard deviation , then the sample coefficient of skewness and coefficient of kurtosis are given by [27–29] The following formula shows the response of the analytical expression of the sample coefficient of skewness: The above study is very important in groundwater study because, to have a clear knowledge of aquifer parameters, several measurements of each parameter must be done, and once these parameters are known, they can be exposed to aleatory uncertainty analysis. In order to point out the possible influence or effect of aperture, in Figures 15 and 16, we present the vertical parallel model with different aperture.
6. Effect of Uncertainty of Aperture on the Advection Dispersion Equation
The consequence of a tracer test of sets of concentration of data is obtained at one or more observation wells or at pumping well at discrete time [30]. The analysis of data starts out from a solution of the analytical or numerical transport equation and determines a set of unknown transport parameters appearing in the solution such that the derivation between measured and calculated concentration values becomes minimal in the least-squares sense [31]. The choice of the solution together with the determined set of parameters constitutes an interpretation of the tracer test data.
The condition of permanent injection is rarely achieved for in situ tracer due to the cost of large amount of tracers and difficulty of monitoring a constant injection flow rate. But it is the first approach to real pollution plumes which are generally detected only after a long period of inflow and subsequent transport under natural aquifer conditions [31]. The corresponding solution is obtained by convolution of the Dirac-input solution. Options for constant dispersivities in flow time are given as follows: is the injection rate of the fluid, is the concentration of injected tracer fluid, is the thickness of the aquifer, is the effective porosity, is the velocity of the fluid, are the longitudinal dispersivities, and is the radioactive decay.
However in the case under investigation, we are not dealing with the radioactive pollution meaning in this case ; also we consider only -direction; then the above equation can be reduced to
The above equation will be called the one-dimensional uniform flow with permanent injection. The one-dimensional solution equation (19) we then used to investigate the effect of uncertainties of the selected injection rate of the tracer fluid is presented earlier in the previous section.
In order to view the numerical results of the above equation, we have made use of the selected 26-injection rate presented earlier as corresponding injection rate from the selected apertures. In addition to this, we chose a typical effective porosity to be 0.05, the longitudinal dispersivities to be 75, the ratio of the borehole to be 0.08 m, the velocity to be 5.1 m/day, and finally the initial concentration to be 100 percent. The graphical representation is depicted in Figure 17 for a fixed distance of 5 m. From Figure 17, it is very important to realize that the size of the aperture plays an important role while dealing with the investigation of the plume. Each size of aperture gives different value of the injection rate and this value of injection gives a different plume; this is depicted in Figure 17.
We present in Figures 18, 19, and 20 some contour plots of some concentrations as function of time and space. This is depicted in Figures 18, 19, and 20.
7. Conclusion
In the recent decade, the idea of fracking has become more and more pronounced in the developed countries. The real motive behind this fracking is the extraction of shale gas which is one of the most useful natural gasses. In our study, we used the mathematical equations describing the relationship between the apertures and the transmissivity, the seepage velocity, and the discharge rate, respectively, to investigate the influence of different apertures in the migration of the pollution through the geological formation called aquifers. To achieve this, we presented a relatively new sampling method, which combines both the Monte Carlo and the Latin hypercube sampling. The method was used to investigate possible risks and uncertainties associated with the process of fracking to extract shale’s gas especially in the case of Karoo system in South Africa. The study reveals that, in the case of the Karoo, the idea of fracking will be successful if and only if a well and the entire fracked reservoir are plugged with, for example, cement, otherwise many aquifers in the Karoo will be polluted.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.