Characterisation of the Hydrogeological Properties of the Ntane Sandstone Aquifer Using Co-Seismic and Post-Seismic Groundwater Level Responses to the Mw 6.5 Moiyabana Earthquake, Central Botswana

Abstract

The 3 April 2017 Mw 6.5 Moiyabana earthquake (Central Botswana) had a significant impact on groundwater levels; a gradual co-seismic increase and a stepwise decline in groundwater levels were observed in response to the earthquake at boreholes MH2 and Z12836, respectively. In this study, we investigated the response of groundwater levels to Earth tides by computing the amplitude and phase shift of the M₂ tidal constituent to estimate the temporal variations of the storativity, transmissivity, and permeability of the Ntane sandstone aquifer (the main aquifer system) prior to and after the earthquake event. The storativity and permeability computed for borehole MH2 showed a decrease in magnitude of 3.17432 × 10⁻⁴ and 1.85 × 10⁻¹³ m² respectively, indicating that strong ground shaking at borehole MH2 might have consolidated the aquifer material, thus resulting in decreased aquifer permeability. The aquifer coefficient of storativity decreased by 2.85 × 10⁻⁴ while permeability was enhanced by 0.047 × 10⁻¹³ m² at borehole Z12836. Enhanced permeability might have resulted from increased/enhanced fracturing of the aquifer, fracture clearing and dynamic shaking.

1. Introduction

Aquifers are responsive to external stresses and periodic loadings such as atmospheric pressure variations, Earth tides and seismic waves [1,2,3,4]. Changes in groundwater levels in confined aquifers occur in response to Earth tides and atmospheric pressure loads [1,5,6,7,8]. Variabilities in well levels caused by earthquakes have been measured and widely studied [4,9,10,11,12]. These variations are sufficiently large to be measured and recorded, and their analysis is very useful in determining the elastic properties and porosity of aquifers [1,13].

Borehole hydrographs are analysed following an earthquake event to identify pre-seismic, co-seismic, and post-seismic groundwater level changes and the observed changes are used to investigate temporal changes in the aquifer’s hydraulic properties [11]. Co-seismic groundwater level responses may vary depending on factors such as the earthquake magnitude, distance from the earthquake epicenter, the magnitude of the seismic energy density, hydrogeological properties of the aquifer system and well construction and design [3,11,14]. Earthquake-induced groundwater level changes thus manifest as abrupt step-like changes (rises or falls) in the near-field, gradual changes in the intermediate field and oscillations in the far-field [4,5,6].

Earthquakes can modify aquifer properties and the observed changes in groundwater level may reflect the equilibration of these properties [7]. The passage of seismic waves through aquifers also contributes to changes in the permeability of the well-aquifer system. When the change in the permeability is not reversed after the seismic waves propagate through the aquifer system, there would be major differences in the groundwater levels before and after an earthquake [8]. The change in transmissivity and storage coefficient induced by earthquakes can cause variations in the phase shift and amplitude of different tidal wave constituents [9]. Thus, the groundwater level tidal response reflects the aquifer properties [8,10,11].

Aquifer parameters, such as specific storage and transmissivity, can be derived from the response of groundwater levels to Earth tides [9,12,13]. The amplitude response of the tides is primarily related to the elastic properties and the porosity of the aquifer [12,14]. However, the phase lag between the tidal dilation of the aquifer and the groundwater-level response in wells is mainly affected by the aquifer’s transmissivity and storage coefficient [15]. Therefore, the computed amplitude response and phase shift tidal response are used to monitor the storativity and transmissivity of an aquifer before and after the earthquake, respectively [2,9].

Mechanisms that account for co-seismic groundwater level changes include (1) enhancement of the aquifer’s permeability, (2) co-seismic liquidation or consolidation of loose sediments and (3) co-seismic static strain and changes in the pore pressure [7,8,16,17]. Analyses of the groundwater level tidal response in wells have confirmed that large earthquakes tend to change the permeability of aquifers locally [8,18,19,20]. Permeability is effectively enhanced through development of new fractures in the aquifer and clearing of mineralised fractures during dynamic shaking of the existing fracture system [8]. However, the increased permeability will eventually decrease over time due to hydrogeological and/or biogeochemical processes that re-clog the fractures [7]. Sustained co-seismic groundwater level changes are attributed to the redistribution of stress and strain due to opening of fractures during earthquakes [3]. Sustained co-seismic groundwater level drop may indicate changes in the poroelastic property of the aquifer system [21].

The 3 April 2017 Mw 6.5 Moiyabana earthquake was the second largest earthquake recorded since monitoring began in Botswana [22,23,24,25]. The earthquake induced distinct co-seismic groundwater level responses in a nearby borehole MH2 and a distant borehole Z12836, which are located 40 km and 167 km from the epicentre, respectively. This earthquake event provided an opportunity to study the post-earthquake hydrogeological responses in the boreholes using time series data of the groundwater level monitoring in the area.

Problem Statement

Research directed at analysing the causative mechanism of the Moiyabana earthquake has received a lot of attention e.g., [22,23,26,27,28,29]. The interest was triggered by the earthquake’s occurrence in an area considered part of a non-seismically active region [22,27,28]. There is, however, a notable lack of research on the effect of this earthquake on the groundwater resources of the nearby and regional aquifer (Ntane Sandstone aquifers). The current research is therefore directed at analysing effects of the earthquake on groundwater resources in the nearby and regional aquifer (Ntane Sandstone aquifer) This is important because large earthquakes can cause significant changes in the hydraulic properties of aquifers and aquitards, with far-reaching implications on the groundwater resources of the affected areas [20]. The Ntane sandstone aquifer covers large areas of central and eastern Botswana and is the main source of water to meet the water demands for domestic and industrial usage across many villages and towns in the area. It is therefore critical to investigate the impact of the earthquake on the hydraulic properties and water resources of the Ntane sandstone aquifer.

The objectives of this study were to analyse the borehole hydrographs for pre-seismic, co-seismic and post-seismic groundwater level responses to the Moiyabana (Mw 6.5) earthquake and to compute the tidal amplitude and phase shift of the M₂ tidal constituent to estimate the aquifer’s storativity, transmissivity and permeability prior and following the earthquake.

2. Materials and Methods

2.1. Study Location

The CIC Energy Wellfield (referred to as the Kudumatse area) is a 150 km² wellfield situated to the southwest of Kudumatse village in-southeast Botswana. The borehole Z12836 (23.510° S and 26.699° E) located in this wellfield is approximately 167 km from the Moiyabana earthquake’s epicentre (Figure 1).

The Gope Well field (referred to as the Gope area) covers an area of approximately 45 km² within the Central “Kalahari Game Reserve (CKGR)”. The wellfield belongs to the Ghaghoo Diamond Mine, formerly known as Gope Mine. The remotely located Ghaghoo Diamond mine is approximately 45 km west of the eastern border of the CKGR. It is within this wellfield that borehole MH2 (22.626° S and 24.762° E) is located approximately 40 km from the Moiyabana earthquake’s epicentre (Figure 1).

The two wellfields are part of the Central Kalahari Basin (CKB) which is predominantly filled by Karoo age sediments that are covered by flat terrain of the younger Aeoalian sands of the Kalahari Group [30,31,32]. The CKB is characterized by a semi-arid to arid climate with characteristic cold, dry winters and hot wet summers [32]. Rainfall is highly variable in the CKB area, both spatially and temporally [33]. Nearly all rainfall occurs as convective thunderstorms during the summer season between September and April, with an average of about 380–530 mm/annum [32].

Figure 1. Map of Botswana showing the Central Kalahari Basin and the major structures associated with it and the location of borehole MH2 and Z12836 (modified from [23,32,34]).

Figure 1. Map of Botswana showing the Central Kalahari Basin and the major structures associated with it and the location of borehole MH2 and Z12836 (modified from [23,32,34]).

2.1.1. Geology

The CKB area is underlain by the Carboniferous to Jurassic age sedimentary succession of the Karoo Supergroup. The Karoo unconformably overlie the Archaean basement and Proterozoic rocks [31,33,34,35]. The Archean basement consists of granitic gneisses of varying compositions, metasediments and the amphibolite to the Limpopo-Shashe Belt [36,37]. The Karoo Supergroup succession is capped by widespread continental flood basalts of the Stormberg Basalts [38]. Younger Tertiary to Quaternary Kalahari sands overlie the CKB rocks and the bedrock geology [31,34,35]. The CKB boundaries have been defined by the E–W-trending Zoetfontein Fault, Limpopo Mobile Belt, Okavango Dyke Swarm to the North and the N–S-trending Kalahari Line to the western edge of the Kaapvaal Craton and the NE–SW-trending Makgadikgadi Line that defines the north-western edge of the Zimbabwe Craton [30,33,39].

Pre-Karoo igneous intrusions that include dolerite dikes and inclined sills have been mapped in the study area [31,32]. The ENE–WSW regional Zoetfontein fault (Figure 1) formed in the Lower Proterozoic age is the most prominent structure. The intrusions emanate from this fault zone along with other conjugate structures due to frequent lateral movement [32,34,40]. Post-Karoo dolerite dykes, sills and diamondiferous kimberlite intruded the Karoo rocks along a WNW–ESE trend [34]. The Zoetfontein Fault was affected by the NW-to-NNW- and NE-trending post-Karoo normal faults [32]. Several faults that trend in the NNW–SSE, ENE–WSW and WNW–ESE orientations transect the Zoetfontein fault, indicating that they post-date it [31,32]. A majority of the faults in the area are younger or the same age as the Stormberg Basalt, as they cut through the entire stratigraphic sequence except for the Kalahari Sand [40,41,42,43]. These normal faults have created a complex set of horst and graben structures [32].

2.1.2. Hydrogeology and Aquifer Parameters

Groundwater in the study area has been limited to the Stormberg basalt, Ntane sandstone and sedimentary sequences of the Ecca Group [30]. The most productive and prominent aquifer in the study area is the Ntane aquifer (comprised of fine-medium-grained sandstone). The Ntane Sandstone aquifer is a dual-porosity aquifer, with the interstitial porosity associated with the poorly cemented Ntane Sandstone Formation and the secondary porosity associated with post-Karoo lineaments that cut through the Karoo lithologies [40].

Lineaments in the Ntane sandstone provide additional pathways for water movement, with the matrix inherently providing matrix storage potential especially where the Ntane sandstone is weathered [40]. The Ntane sandstone aquifer has high transmissivity values where fracturing is well developed and cementation is less, hence the Ntane sandstone’s high porosity and permeability are attributed largely to intense fracturing [33]. The Ntane sandstone aquifer is confined at the top by the basalt cap and below by the low permeability mudstone and siltstone of the Mosolotsane Formation (Figure 2a,b). Where the Stormberg basalt cap is missing or completely eroded, the Ntane sandstone is hydraulically connected with the overlying Kalahari Sands and thus constitutes an unconfined aquifer [40]. The regional groundwater rest levels are approximately 100 m and 34 m in the Gope and Kudumatse areas respectively, indicating confined aquifer conditions [41,42,43].

The hydraulic conductivity (K), transmissivity (T) and Storativity (S) of the Ntane aquifer were calculated using the Cooper-Jacob (1946) technique under a constant-rate pumping test [44]. The results are as follows; for borehole MH2, K = 0.098 m/day, T = 6.822 m²/day and S = 0.0001, and for borehole Z12836, K = 0.488 m/day T = 33.696 m²/day and S = 0.0002 [41,42,43]. These storativity values are indicative of confined aquifer conditions. The Kudumatse area is considered a high transmissivity zone, with T ranging from 500–3300 m²/day, and this is attributable to heavy fracturing and the low cementation of the sandstone [42].

2.1.3. Borehole Construction

Neither of the boreholes were screened, and thus the steel casing installed in the boreholes was perforated at the aquifer interval to allow the ingress of groundwater into the borehole. The steel casing for borehole Z12836 had a diameter of 254 mm and 203 mm from intervals 0–17 mbgl and 1 to 187 mbgl, respectively. The wellbore is open at 104–173 mbgl at Ntane sandstone, which is the principal aquifer. The perforated casing length is 69 m. It should be noted that the Ntane sandstone is friable towards the top of the stratigraphic succession, and it becomes moderately competent towards the bottom.

The casing diameter for borehole MH2 for example was 203 mm. The wellbore is open at 394–425 mbgl at Ntane sandstone aquifer, which is the main water bearing rock (

Table 1. Summary of the hydrogeological conditions at boreholes MH2 and Z12836.

Borehole ID	Borehole Depth (mbgl)	Depth to Base (mbgl)				Water Strike (mbgl)	Rest Water Level (mbgl)
		Kalahari Beds	Stormberg Basalt	Ntane Sandstone Fm	Mosolotsane Fm
MH2	501	71.8	394	425	>501	211; 388	99.362
Z12836	187	26	104	173	>187	58; 104; 112	34.19

). Therefore, the perforated casing length is 31 m and the borehole drilling diameter was 254 mm.

2.2. Methodology

2.2.1. Sampling

Groundwater level was recorded at a temporal resolution of 6 h at boreholes Z12836 and MH2. Barometric pressure was measured for borehole BH7435 at a temporal resolution of 12 h using the Baro-Diver data logger, which recorded the semi-diurnal atmospheric pressure variations. Borehole BH 7435 is situated 25 km from borehole Z12836. The groundwater level and barometric pressure measurements were recorded during the period of February 2016 to March 2018 at borehole Z12836, while at borehole MH2, the groundwater level was recorded from January 2016 to December 2018.

2.2.2. Borehole Hydrograph Analysis

Borehole hydrographs were analysed to identify the pre-, co- and post-seismic groundwater level responses due to the Moiyabana earthquake. Earthquake-induced groundwater level changes can be detected by analysing anomalous signals of groundwater level time series information from a monitoring well, in order to detect earthquake-induced groundwater level changes [1]. Co-seismic groundwater level changes manifest themselves as step-like changes, gradual changes and oscillations [5,45]. Post-seismic water level adjustment or recovery is termed ‘post-seismic groundwater recession’ [4,46].

2.2.3. Tidal Analysis

After all the data were resampled into hourly measurements, we applied a least-squares method for tidal analyses using the Tsoft software program, Version 2.2 [9,47]. Spectral analysis was performed on the MH2 and Z12836 groundwater levels, and the barometric pressure time series data were used to identify periodic fluctuations in the data and to extract the dominant tidal components using the Tsoft software package [47]. Tsoft uses a Fast Fourier Transformation (FFT) algorithm to calculate time series data spectra [48,49].

A low-pass filter and a high-pass Butterworth filter with a cut-off frequency equal to 3.0 cpd and 0.8 cpd (cycles per day) respectively were applied to the groundwater level data to enhance the signal-to noise ratio of the main tidal constituents [10,50]. In this study, barometric loading was extracted from the groundwater level data using the Bayesian Tidal Analysis Program (Baytap-G) software, Version 1.9.2 developed in [51]. This was achieved using the non-harmonic method based on multiple linear regression [52]. Borehole MH2 had no barometric pressure record; therefore no barometric pressure correction was performed on the groundwater levels data. The impact of the lack of barometric pressure record on the M₂ tidal component of the groundwater level for borehole MH2 is likely to have minimal impact on the M₂ component of the groundwater level [9,21]. Moreover, since borehole MH2 and Z12836 are far from the coastline, ocean tide corrections were also deemed not very necessary [21].

The main tidal wave components that have a significant effect on groundwater level variations are O_1, M₂, N_2, K₁ and S₂ waves,

Table 2. Major harmonic components of the tidal potential (adapted with permission from Munk& MacDonald [53], 1960, Munk & MacDonald).

Tidal Component	Description	Period (Day)	Frequency, cpd (Cycles per Day)
O1	Principal lunar	1.0758	0.9295
K1	Lunar–solar	1.3721	1.0029
N2	Lunar elliptic	0.5275	1.8957
M2	Principal lunar	0.5175	1.9324
S2	Principal solar	0.5000	2.0000

[2,6,53]. Tidal waves greater in amplitude of the five major tidal components, i.e., M₂ and O₁, are usually preferred to evaluate the groundwater-level response to earth tides because they have a greater signal-to-noise ratio [6,9,49]. However, in this study, only the lunar semi-diurnal (M₂) wave constituent was used to estimate the aquifer hydraulic properties. This is because it is more stable, exhibits a large amplitude and the smallest root-mean-square error, and is less affected by external influences from solar radiation/thermal effect and diurnal barometric pressure effects than the O₁ tidal component [2,9,10,11,21,46,49].

A window size of 31 days was chosen for this study to separate the time domain tidal analysis for the semidiurnal M₂ and the diurnal S₂ tidal components and to avoid spectral leakage [10,54,55]. Thus, the response of groundwater level to the volumetric strain associated with the M₂ tidal constituent was obtained every 31 days using the Baytap-G software. The pre-earthquake period was taken as 31 days before the earthquake and 31 days after the earthquake for the post-earthquake period, allowing for a comparison to be made, and also to be determine the impact of the Moiyabana earthquake on the hydraulic properties of the aquifer. In order to estimate earthquake-induced changes in the hydraulic properties of the Ntane sandstone aquifer system, the groundwater level response to earth tides and atmospheric pressure before and after the Moiyabana earthquake’s occurrence were used.

2.3. Mechanisms for Classifying/Determining the Hydrogeological Properties

2.3.1. Estimating the Aquifer Storage Coefficient

The variation in groundwater level produced by aquifer dilatation as a result of Earth tide is a function of the aquifer’s specific storage and can be calculated by measuring fluctuations in groundwater level (dh) and tidal component amplitude [9]. The storage coefficient $\mathrm{S}$ is the product of the aquifer specific storage $\mathrm{S}\mathrm{s}$ and the thickness of the aquifer, $\mathrm{d}$. The specific storage $\mathrm{S}\mathrm{s}$ therefore can be calculated from the groundwater level tidal-response amplitude [6,56,57] using Equation (1). The minus sign in Equation (1) indicates that the tide-generating potential increases when the head in the aquifer decreases [57].

$$\mathrm{S}\mathrm{s}=\left\{\left[\left(\frac{1-2\mathrm{v}}{1-\mathrm{v}}\right)\left(\frac{2\overline{\mathrm{h}}-6\overline{\mathrm{l}}}{\mathrm{a}}\right)\right]\frac{\mathrm{d}{\mathrm{W}}_2}{\mathrm{g}}\right\}\left({-\mathrm{d}\mathrm{h}}^{-1}\right)$$

(1)

where $\mathrm{v}$ is the Poisson ratio, $\overline{\mathrm{h}}$ is the love number of the surface of the earth (0.60), $\overline{\mathrm{l}}$ is the love number at surface of the earth (0.07) [53], $\mathrm{a}$ is the mean earth radius $6.371\times {10}^8 \mathrm{c}\mathrm{m}$ and ${\mathrm{W}}_2$ is the theoretical tidal potential ${\mathrm{W}}_2\approx {\mathrm{A}}_2\left(\mathsf{\theta },\mathsf{\delta }\right)=\mathrm{g}{\mathrm{K}}_{\mathrm{m}}\mathrm{b}\mathrm{f}\left(\mathsf{\theta }\right)$ [57]. The terms enclosed in the brackets will always be constant for general aquifer materials i.e., 7.88 × 10⁻⁹, hence:

$$\mathrm{S}\mathrm{s}=7.880\times {10}^{-9}·\frac{{\mathrm{A}}_{2(\mathsf{\tau },\mathsf{\theta })}}{{\mathrm{A}}_{\mathrm{h}}\left(\mathsf{\tau }\right)}$$

(2)

where is ${\mathrm{A}}_{\mathrm{h}}\left(\mathsf{\tau }\right)$ the fluctuation amplitude induced by the different tidal period components; ${\mathrm{A}}_{2(\mathsf{\tau },\mathsf{\theta })}$ is the fluctuation amplitude of the second-step tide-generating potential; ${\mathrm{K}}_{\mathrm{m}}$ is the general lunar coefficient, typically taking a value of ~0.537 m [53]; $\mathrm{b}$ is a constant related to the period of the various tidal wave component or the amplitude factor for various tidal components; and $\mathrm{f}\left(\mathsf{\theta }\right)$ is a function of latitude $\mathsf{\theta }$ or the latitude dependence factor. For the semi-diurnal M₂ wave, $\mathrm{f}\left(\mathsf{\theta }\right)=0.5{\mathrm{c}\mathrm{o}\mathrm{s}}^2\left(\mathsf{\theta }\right)$ and $\mathrm{b}\approx 0.908$ [53,57].

After obtaining the specific storage $\mathrm{S}\mathrm{s}$ using Equation (2), the storage coefficient $\mathrm{S}$ can be calculated using $\mathrm{S}=\mathrm{S}\mathrm{s}\times \mathrm{d}$. d is the thickness of the aquifer.

2.3.2. Estimation of the Aquifer Transmissivity Using Phase Shift

The Horizontal Flow Model

The Hsieh solution [9] was used to calculate the transmissivity of an aquifer when the phase shift of the groundwater level relative to the tidal strain is negative i.e., ranging from −90° to −1° [58]. The Hsieh solution is used for the tidally induced radial flow of groundwater to a well (assuming a horizontal, laterally extensive, confined, isotropic and homogeneous aquifer that is fully penetrated by a well) to estimate the aquifer transmissivity from the phase shift associated with the M₂ tidal constituent. The Hsieh solution can be approximated for a realistic well geometry and aquifer properties as

$${\mathrm{s}}_{\mathrm{w}}=\frac{{\mathrm{Q}}_0}{2\mathsf{\pi }\mathrm{T}}\left\{\left[\mathsf{\varphi }\mathrm{K}\mathrm{e}\mathrm{r}\left(\mathsf{\alpha }\mathsf{\omega }\right)-\mathsf{\psi }\mathrm{K}\mathrm{e}\mathrm{i}\left(\mathsf{\alpha }\mathsf{\omega }\right)\right]+\mathrm{i}\left[\mathsf{\psi }\mathrm{K}\mathrm{e}\mathrm{r}\left(\mathsf{\alpha }\mathsf{\omega }\right)+\mathsf{\varphi }\mathrm{K}\mathrm{e}\mathrm{i}\left(\mathsf{\alpha }\mathsf{\omega }\right)\right]\right\}$$

(3)

With

$$\mathsf{\varphi }=\frac{-\left[{\mathrm{K}\mathrm{e}\mathrm{r}}_1\left(\mathsf{\alpha }\mathsf{\omega }\right)+{\mathrm{K}\mathrm{e}\mathrm{i}}_1\left(\mathsf{\alpha }\mathsf{\omega }\right)\right]}{2^{1/2}\mathsf{\alpha }\mathsf{\omega }\left[{\mathrm{K}\mathrm{e}\mathrm{r}}_1^2\left(\mathsf{\alpha }\mathsf{\omega }\right)+{\mathrm{K}\mathrm{e}\mathrm{i}}_1^2\left(\mathsf{\alpha }\mathsf{\omega }\right)\right]}$$

(4)

$$\mathsf{\psi }=\frac{-\left[{\mathrm{K}\mathrm{e}\mathrm{r}}_1\left(\mathsf{\alpha }\mathsf{\omega }\right)+{\mathrm{K}\mathrm{e}\mathrm{i}}_1\left(\mathsf{\alpha }\mathsf{\omega }\right)\right]}{2^{1/2}\mathsf{\alpha }\mathsf{\omega }\left[{\mathrm{K}\mathrm{e}\mathrm{r}}_1^2\left(\mathsf{\alpha }\mathsf{\omega }\right)+{\mathrm{K}\mathrm{e}\mathrm{i}}_1^2\left(\mathsf{\alpha }\mathsf{\omega }\right)\right]}$$

(5)

$$\mathsf{\alpha }\mathsf{\omega }={\left(\mathsf{\omega }\mathrm{S}/\mathrm{T}\right)}^{1/2}·{\mathrm{r}}_{\mathrm{W}}$$

(6)

$$\mathsf{\omega }=2\mathsf{\pi }/\mathsf{\tau }$$

(7)

${\mathrm{s}}_{\mathrm{w}}$ is the drawdown at the well, which is related to the water level in the well (x) and the pressure head in the aquifer (h) by h—${\mathrm{s}}_{\mathrm{w}}$ = x, Q₀ is the discharge of the aquifer at the well, $\mathrm{K}\mathrm{e}\mathrm{r}\left(\mathsf{\alpha }\mathsf{\omega }\right)$ and $\mathrm{K}\mathrm{e}\mathrm{i}\left(\mathsf{\alpha }\mathsf{\omega }\right)$ are the Kelvin functions of order zero and ${\mathrm{K}\mathrm{e}\mathrm{r}}_1\left(\alpha \omega \right)$ and ${\mathrm{K}\mathrm{e}\mathrm{i}}_1\left(\mathsf{\alpha }\mathsf{\omega }\right)$ are the Kelvin functions of order one. $\mathsf{\omega }$ is the tidal fluctuation frequency of the tidal wave constituent, taking the value of $\mathsf{\tau }$ to be 0.5175/d for the M₂ wave. ${\mathrm{r}}_{\mathrm{w}}$ is the inner radius of the well casing. T and S are the transmissivity and storativity of the aquifer, respectively.

The relationship between aquifer properties and the amplitude ratio, which is the ratio of the amplitude of tide of the water level in a well to the amplitude of the fluctuating pressure head in the elastic aquifer responding to the tidal stress, can be obtained [9]. The amplitude response A and the phase shift $\mathsf{\eta }$ can be defined using Equations (8) and (9) respectively.

The amplitude response $\mathrm{A}$ is:

$$\mathrm{A}={({\mathrm{E}}^2+{\mathrm{F}}^2)}^{-1/2}$$

(8)

While the phase shift $\mathsf{\eta }$ is:

$$\mathsf{\eta }=-a\mathrm{r}c{tan}^{-1}(F/E)$$

(9)

where E and F are defined as

$$\mathrm{E}=1-\frac{{{\mathsf{\omega }\mathrm{r}}_{\mathrm{c}}}^2}{2\mathrm{T}}\left[\mathsf{\psi }\mathrm{K}\mathrm{e}\mathrm{r}\left(\propto \mathsf{\omega }\right)+\mathsf{\varphi }\mathrm{K}\mathrm{e}\mathrm{i}\left(\propto \mathsf{\omega }\right)\right]$$

(10)

$$\mathrm{F}=\frac{{{\mathsf{\omega }\mathrm{r}}_{\mathrm{c}}}^2}{2\mathrm{T}}\left[\mathsf{\varphi }\mathrm{K}\mathrm{e}\mathrm{r}\left(\propto \mathsf{\omega }\right)-\mathsf{\psi }\mathrm{K}\mathrm{e}\mathrm{i}\left(\propto \mathsf{\omega }\right)\right]$$

(11)

${\mathrm{r}}_c$ is the inner radius of screened or open portion of the well. For Earth tide analysis, realistic values of ${\mathrm{r}}_{\mathrm{w}}$, T and S, and the value of $\propto \omega$ computed by Equation (6) will usually be small (<0.1). In this case, both ${\mathrm{K}\mathrm{e}\mathrm{r}}_1\left(\mathsf{\alpha }\mathsf{\omega }\right)$ and ${\mathrm{K}\mathrm{e}i}_1\left(\mathsf{\alpha }\mathsf{\omega }\right)$ can be approximated by $-1/\left({2^{1/2}}_{{\propto }_{\mathsf{\omega }}}\right)$. This leads to $\mathsf{\varphi }\approx 0$ and $\mathsf{\psi }\approx 0$ and E and F can thus be approximated according to Equations (12) and (13) respectively as:

$$\mathrm{E}\approx 1-(\mathsf{\omega }{\mathrm{r}}_{\mathrm{c}}^2/2\mathrm{T})\mathrm{K}\mathrm{e}\mathrm{i}\left({\mathsf{\alpha }}_{\mathrm{w}}\right)$$

(12)

$$\mathrm{F}\approx (\mathsf{\omega }{\mathrm{r}}_{\mathrm{c}}^2/2\mathrm{T})\mathrm{K}\mathrm{e}\mathrm{r}\left({\mathsf{\alpha }}_{\mathrm{w}}\right)$$

(13)

Using ${\mathrm{r}}_c$ of 0.1015 m for both boreholes, τ for the M₂ tidal constituent is 12.421 h (44,715.6 secs); the transmissivity was calculated using the observed phase shift taking the well-aquifer parameter $\mathrm{S}{\mathrm{r}}_{\mathrm{w}}^2/{\mathrm{r}}_{\mathrm{c}}^2=1\times {10}^{-4}$ or negative phase shift as illustrated in [9].

Vertical Flow Model

The vertical flow model is used to estimate the transmissivity of an aquifer when the computed phase shift is positive. The vertical flow model ignores horizontal flow and predicts the phase shifts as positive (0–45°) between the groundwater level and the tidal strain [2,59]. The vertical flow model cannot cause negative phase lag beyond −1°. It is satisfied by undrained boundary conditions at an infinite depth and drained at the water table. Thus, the relationship between the groundwater level response inside a wellbore and the tidal fluctuations under these conditions can be estimated using Equations (14) and (15) [60].

The amplitude response $\mathrm{A}$ is defined as;

$$\mathrm{A}=\left|\frac{{\mathrm{x}}_0}{{\mathsf{\epsilon }}_0}\right|=\frac{1}{\mathrm{S}\mathrm{s}}{\left[1-2\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{\mathrm{z}}{\mathsf{\delta }}\right)\cos \left(\frac{\mathrm{z}}{\mathsf{\delta }}\right)+\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{2\mathrm{z}}{\mathsf{\delta }}\right)\right]}^{-1/2}$$

(14)

The phase shift $\mathsf{\eta }$ is defined as

$$\mathsf{\eta }=\mathrm{a}\mathrm{r}\mathrm{g}\left(\frac{{\mathrm{x}}_0}{{\mathsf{\epsilon }}_0}\right)={\tan }^{-1}\left\{\frac{\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{\mathrm{z}}{\mathsf{\delta }}\right)\sin \left(\frac{\mathrm{z}}{\mathsf{\delta }}\right)}{1-\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{\mathrm{z}}{\mathsf{\delta }}\right)\cos \left(\frac{\mathrm{z}}{\mathsf{\delta }}\right)}\right\}$$

(15)

The empirical variable is $\mathsf{\delta }=\sqrt{2\mathrm{D}/\mathsf{\omega }}$, where $\mathrm{D}$ is the hydraulic diffusivity which is defined as the ratio of Transmissivity to Storativity and $\mathrm{z}$ is the depth from the surface to the bottom of the borehole. Thus, the transmissivity of the aquifer can be estimated using the measured phase shift by Equation (16).

2.3.3. Estimation of the Aquifer’s Permeability

The transmissivity $\mathrm{T}$ of an aquifer is estimated by following the standard method of fitting the tidal data for phase shift $\mathsf{\eta }$ using Equations (11) and (15), and using the measured values of r_c and r_w and the approximated value of S [8]. Then the estimated transmissivity is used to estimate the permeability of the aquifer. The relationship between $\mathrm{T}$ and permeability, $\mathrm{k}$ can be obtained using Equation (16) as follows:

$$\mathrm{k}=\frac{\mathsf{\mu }\mathrm{T}}{\mathrm{b}\mathsf{\rho }\mathrm{g}}$$

(16)

b denotes the thickness of the aquifer. $\mathsf{\rho }$ is the density, $\mathsf{\mu }$ is the coefficient of dynamic viscosity of water, and g the gravitational acceleration and the typical values are 10⁻³ Pa·s, 10³ kg/m³ and 9.81 m/s², respectively [61].

3. Results

3.1. Borehole Hydrograph Analysis

A gradual upward co-seismic groundwater level displacement of 1.565 m was observed in the hydrograph of borehole MH2 (Figure 3). This was followed by a post-seismic groundwater level recession. A full recovery to the pre-seismic groundwater levels was achieved after fifteen months (i.e., July 2018). The groundwater level recovery is also related to the permeability and thickness of the aquifer, with a thinner the aquifer, yielding a faster the rate of recovery of the co-seismic groundwater level change [3]. The recovery of the groundwater level to the pre-seismic groundwater level may show that re-equilibration of the aquifer properties occurred, which was changed by the earthquake’s occurrence. Borehole MH2 is a dewatering borehole, as the mine is under care and maintenance. The anomalous rise at the end of the hydrograph might be indicating non-pumping periods, and therefore, recovery towards natural levels.

A 0.2201 m co-seismic step-like fall in the groundwater level was observed in the hydrograph of borehole Z12836. Recovery to the pre-seismic groundwater level was not achieved in the period considered in this study (up to April 2018) (Figure 4). The co-seismic groundwater level was sustained in borehole Z12836, which might be an indication of a permanent effect on the aquifer. As there was lack of groundwater level recovery observed at borehole Z12836 after the Mw 6.5 Moiyabana earthquake, this may be suggesting the occurrence of groundwater discharge, possibly through fractures generated by the earthquake [3,62]. As seismic shaking causes fracturing of the aquifer, the groundwater outflow through the fractures would further reduce pore pressure. Consequently, the aquifer can no longer return to the original state [3]. The earthquake triggered fractures might have created conduits through the confining layer to dissipate the pore pressure of the aquifer [3]. The confinement of the aquifer may also be breached by the fracturing due to seismic shaking during an earthquake. The dissipation of excess pressure in the confined aquifer generated by earthquakes is a slow process due to the retardation of the less permeable confining layer [3]. Therefore, groundwater level changes in confined aquifers are likely to be sustained for a long time.

Seismic energy density in excess of 10⁻³ J/m³ causes sustained groundwater level changes in wells. Seismic energy density is the maximum seismic energy required in a unit volume of rock or sediment to do work [4]. Thus, the seismic energy density $\mathrm{e}$ (J/m³) is calculated according to the empirical relationship between earthquake magnitude $\left(\mathrm{M}\mathrm{w}\right)$ and earthquake epicentral distance $\mathrm{r}$ (km) [4] as follows:

$$\log \mathrm{r}=0.48\mathrm{M}\mathrm{w}-0.33\log \mathrm{e}-1.4$$

(17)

The magnitude of the Moiyabana earthquake was Mw 6.5 and

Table 3. The computed seismic energy density experienced at the boreholes.

Borehole ID	Epicentral Distance (km)	Seismic Energy Density (J/m3)
MH2	40	2.277
Z12836	167	0.0299

indicates the epicentral distance used for the seismic energy density computations. The calculated seismic energy density values using Equation (17) for borehole MH2 and Z12836 are 2.277 J/m³ and 0.0299 J/m³ respectively (Table 3), which are within the range of seismic energy density capable of triggering sustained groundwater level changes induced by earthquakes; these values are also consistent with previous studies using global datasets [4,63]. Similarly, most sustained groundwater level changes are bounded by $\mathrm{e} \times$ 10⁻³ J/m³ [4].

3.2. Spectra Analysis

The groundwater level spectra of the MH2 borehole show that the amplitude of the S₂ tidal component is greater than the M₂ tidal component (Figure 5). Therefore, because of this tidal component’s low signal-to-noise ratio, the interpretation of the tidal response of the M₂ tidal component at borehole MH2 requires more caution [55].

The computed groundwater level spectra for borehole Z12836 show O₁, K₁, N₂, M₂ and S₂ tidal constituents and indicate excellent tidal responses (Figure 6). The spectrum of the groundwater level indicates that fluctuations in the groundwater level in the boreholes occur at the same frequency as fluctuations in the pressure head in the aquifer induced by the Earth tides [49]. The N₂, component is small, and shows at a slightly lower frequency than the M₂ component (Figure 6). Groundwater level data associated with inland wells show the presence of the lunar harmonics O₁ and M₂ representing the effect of Earth tides, as the influence of ocean tides should be negligible and lunar harmonics do not have periodic oscillations in atmospheric pressure [57]. At a frequency of 1.9324 cpd corresponding to the M₂ tidal component, only Earth tide influences are present [48].

Clear atmospheric tidal harmonics at 1 cpd (S₁) constituent are shown by the barometric pressure spectra, which occur at the same harmonic frequency as the K₁ observed in the barometric spectra in Figure 7 [9,10,48]. The atmospheric harmonic S₁ can therefore interfere with the water level tidal reaction at this frequency [58]. Therefore, to estimate the hydrogeological properties, this study uses the stable M₂ constituent, which is present in the groundwater level spectrum of the two boreholes (MH2 and Z12836), as it is less affected by barometric pressure loading contamination [55].

3.3. Tidal Analysis

Seismic waves often change the permeability and storage properties of an aquifer system. Therefore, these properties can be estimated from the amplitude responses and phase shifts of tidal components of the groundwater level [8,15,64]. Thus, temporal variations in the tidal constants might reflect changes in aquifer hydraulic properties. In order to investigate potential earthquake-related temporal variation in hydraulic parameters, we measured the variability in amplitude and phase changes in groundwater level fluctuation relative to the M₂ tidal constituent’s volumetric strain. Hsieh et al. [9] indicated that the earthquake-induced changes in transmissivity and storage coefficient can cause phase shift and amplitude variations of the tidal constituents.

The tidal amplitude and phase variation in the groundwater level to Earth tides for borehole MH2 are illustrated in Figure 8. The tidal amplitude ranged from 1.5 mm to 3.5 mm for borehole MH2 (Figure 8a). The observed amplitude response is constant without major variations, indicating that the aquifer’s storage properties are constant and display no substantial changes [58]. The amplitude did not show major deviations except for April 2017, when the Moiyabana earthquake took place. The Moiyabana earthquake greatly enhanced the amplitude response of the groundwater level for borehole MH2. The tidal constants remained stable after the earthquake’s occurrence. However, the changes in amplitude can be attributed to changes in storativity of the aquifer [8,9]. The amplitude variation of the M₂ was very small, which might be an indication that the storativity of the aquifer did not change significantly except after the Moiyabana earthquake.

The phase lag obtained from the tidal response of the groundwater level ranged from 75° to 117° and for borehole MH2 (Figure 8b) over the period considered (January 2016–December 2018) in this study with no major deviations. The positive and negative phase shifts indicate the groundwater level advancing or lagging behind the tidal strain, respectively [10]. The computed t phase shift of the groundwater level for the M₂ tidal constituent in boreholes MH2 were both significantly changed by the Moiyabana earthquake (Figure 8b). The phase shift change from 10.592° to −3.434° after the earthquake took place, as observed in borehole MH2 (Figure 8b), may be indicating a decrease in the permeability of the aquifer.

For borehole Z12836, the phase shift changed from −2.933° to positive 12.363° after the earthquake’s occurrence, suggesting that the permeability of the well-aquifer system may have been enhanced by the seismic waves (Figure 9). Phase shift changes have been documented to reflect changes in the aquifer permeability due to unclogging of pre-existing cracks by seismic waves [8]. In a confined aquifer, small phase lags are due to high permeability, while large phase lags are due to low permeability [11]. Thus, an increase in phase shift implies an increase in permeability [8,9]. However, the changes in amplitude can be attributed to changes in storativity of the aquifer [8,9]. After the Moiyabana earthquake, the amplitude ratio and phase shift remained constant, which may suggest that the hydraulic property changes caused by the earthquake at borehole Z12836 had not recovered [8,65] by the end of the analysed time span in this study.

These results are consistent with previous studies [2,6,8], which demonstrate that far-field earthquakes can enhance the aquifer permeability. The results of [2] showed that both the amplitude of the tidal response and the phase shift of the M₂ constituent increased, and that the post-earthquake phase shift changed from negative to positive and an increased in permeability was observed in the aquifer. In [8], characteristic M₂ wave co-seismic and post-seismic phase shifts observed in the wells in California suggested an enhancement of the aquifer’s permeability as a potential cause for the phase shifts. Therefore, permeability enhancement following an earthquake event may be related to propagation of seismic waves.

Previous studies indicated that by analysing the tidal phase and amplitude ratio of the tidal constituents, the effect of enhanced horizontal permeability and increased vertical permeability could be distinguished [2,9,54,58]. A substantial change in the phase shift, along with in decreased amplitude ratio, suggests an enhancement of both the horizontal and vertical permeability [59]. However, the results of the tidal analysis of the Z12836 borehole showed a substantial increase in the phase shift (−2.933° to 12.363°) but a fairly small increase in the amplitude ratio (0.07614 to 0.10879) following the Moiyabana earthquake. A significant increase in the phase shift and small changes in the amplitude ratio can also imply a bidirectional increase in permeability [2,58].

The tidal analysis summary and the estimated hydraulic properties of the Ntane sandstone aquifer at boreholes MH2 and Z12836 before and after the earthquake are shown in

Table 4. Summary of the M2 wave tidal and well-aquifer system hydraulic parameters for borehole. MH2 and Z12836 before and after the earthquake occurrence.

BH ID	Date	Tidal Amplitude (mm)	Phase Shift (°)	Specific Storage Ss (10−6 m−1)	Storativity (Dimensionless) S (10−4)	Transmissivity T (10−5 m2/s)	Permeability k (10−13 m2)
MH2	Mar 17	1.592	10.592	20.172	6.25332	9.7766	3.215
	Apr 17	3.209	−3.437	9.933	3.079	9.2157	3.03
	Δ	1.617	−14.029	−10.521	−3.17432	−0.5609	−0.185
Z12836	Mar 17	2.46	−2.933	12.883	8.89	8.0638	1.191
	Apr 17	3.514	12.363	8.752	6.04	8.3799	1.238
	Δ	1.054	15.259	−4.131	−2.85	0.3161	0.047

. The temporal variations in transmissivity are interpreted as variations in permeability for the duration considered in this study as the aquifer thickness and fluid properties are unlikely to differ. The storativity of borehole MH2 decreased with a magnitude of 3.17432 × 10⁻⁴ after the earthquake. The transmissivity and permeability of borehole MH2 showed a decrease with a magnitude of 0.5609 × 10⁻⁵ m/s and 0.185 m²/day, respectively. The storativity decreased with a magnitude of 2.85 × 10⁻⁴ while the permeability was enhanced by 0.047 m² after the Moiyabana earthquake occurrence for borehole Z12836.

4. Discussion

The storativity values (Table 4) obtained by the method of tidal analysis are consistent with the values obtained from the pumping test. A constant discharge rate aquifer test analysis resulted in a storativity of 1 × 10⁻⁴ and 2 × 10⁻⁴ for boreholes MH2 and Z12836, respectively, using the Cooper-Jacob analysis method. These values are within reasonable ranges of storativity for confined aquifers i.e., 0.5 × 10⁻⁵ to 5 × 10⁻³ [61,66]. Obtaining specific storage estimates within the same order of magnitude as the tidal analysis approach is satisfactory and the findings are consistent with the pumping test results.

The transmissivity estimate for borehole MH2 from the pumping test and tidal analysis were 6.822 m²/day and 8.44 m²/ day, respectively. This indicates that the transmissivity estimate from the tidal analysis method is within an acceptable range and is consistent with the constant rate test result, despite minor deviations. Although the scale of the region sampled by the tidal response is small, the pumping test and the tidal response yield a similar order of magnitude of permeability values. This might be an indication that within the volume investigated by tides, the wells are efficiently interconnected with the fracture network [10].

For borehole Z12836 the constant rate test yielded a transmissivity and permeability of 3.9 × 10⁻⁴ m²/s and 5.899 × 10⁻¹³ m², respectively. The distinction in the obtained estimates of the transmissivity and permeability using the tidal analysis and pumping test could be that pumping test results are affected by the average properties inside the effective radius of the borehole. Therefore, presence of secondary fissures could result in higher hydraulic conductivity values estimates from pump test data analyses [67]. However, for a well-aquifer system that responds to tidal strain, the wells may not efficiently interconnect with the fracture network, which leads to lower transmissivity estimates. Thus, the estimated tidal analysis transmissivity values of borehole Z12836 are slightly lower than the pumping test values by an order of a magnitude.

The estimated permeability values from the tidal analysis are slightly lower than the permeability values from the pumping test estimates due to the influence of a larger skin effect in the tidal response method [58]. The higher initial displacement of the well water level during a pumping test can disrupt the skin effect; thus, a lower skin effect during the pumping test would result in a higher transmissivity value [10,58,68]. Although there is a slight difference in permeability estimates from the two methods, the difference is within a reasonable range. The target aquifer is often not fully confined (semi-confined) therefore it may be characterized by vertical leakage of the confining aquitard, which appears to result in a smaller amplitude ratio and a larger phase advance or a smaller phase lag in the aquifer below the aquitard [68]. An aquitard breach during seismic wave propagation can also cause vertical leakage. The pumping test and the tidal response yield similar permeability values, despite the smaller scale of the region sampled by the tidal response. This indicates that the wells are efficiently interconnected to the fracture network within the volume investigated using Earth tide response [10].

The sampling scale of a pumping test is the distance between the monitoring well and the pumping well [10]. Tidal responses are susceptible to scales ranging from one to tens of meters that are difficult to achieve with traditional methods such as pumping and slug tests [13]. The tidal process averages the properties of aquifers extending an effective radius of the order of $\sqrt{\mathrm{D}\mathsf{\tau }}$ from the well [8,58], where D (the transmissivity to storativity ratio) is the hydraulic diffusivity and $\mathsf{\tau }$ is the tidal period. The approximate values evaluated around borehole MH2 using the tidal responses are in the range from 63 to 78 m and from 58 to 65 m for boreholes MH2 and Z12836, respectively, indicating that the scales of the values based on tidal response are similar to those of the pumping tests. However, larger-scale tests such as pumping tests in fractured-rock aquifers can capture more interconnected zones of fractures and preferential flow paths, which translate into higher transmissivity values [69].

Ground shaking caused by an earthquake may lead to the re-arrangement of solid grains, thus compacting them [70]. As borehole MH2 is located in the intermediate field, the strong ground shaking experienced might have contributed to the consolidation of the aquifer material and a resultant decrease in permeability. Volumetric changes in the rock mass in the near field cause pore pressure changes, thereby causing abrupt groundwater level increase in the nearby wells [71]. Increases in the groundwater level may also result from undrained dilatation or consolidation of sediments caused by strong ground shaking near the earthquake hypocentre. Decreases in groundwater level may be a sign of dilatancy related to fracture formation near the rupture fault [72].

The enhanced aquifer permeability at borehole Z12836 by the Moiyabana earthquake might be due to fracturing of the aquifer material [2,4,8,16,71]. We can infer that fractures which were blocked by colloids before the earthquake were unclogged as there was sufficient seismic energy for groundwater flow which drove colloids away from the fracture. This might lead to an increase in permeability after the seismic wave propagation [8,16]. In addition, fracture clearing and increased permeability due to the consolidation of loose sediments caused by the earthquake-induced dynamic stress have been widely used to explain most far-field water-level changes [2,8]. Dynamic shaking can also enhance effective permeability, especially that of fractured systems [8].

The groundwater levels of borehole Z12836 could not recover to the pre-earthquake water levels (after the Mw 6.5 Moiyabana earthquake); this suggests an increased groundwater discharge possibly through fractures generated or enhanced by the earthquake [3,6]. Seismic shaking is known to cause fracturing of the aquifer, and thus increasing groundwater outflow through the fractures would reduce pore pressure. Consequently, the aquifer in the long run cannot return to the original state [3].

5. Conclusions

The Mw 6.5 Moiyabana earthquake significantly affected the groundwater levels in borehole MH2 and Z12836 and subsequently the hydraulic properties of the Ntane sandstone aquifer. The storativity of the Ntane sandstone aquifer around borehole MH2 and Z12836 decreased after the Mw 6.5 Moiyabana earthquake. Following the earthquake’s occurrence, the permeability was decreased at borehole MH2, while at borehole Z12836, the permeability of the aquifer was enhanced. Changes in aquifer permeability may have a major effect on groundwater flow, and therefore on groundwater sources and the potential for solute transport may increase; thus, they can be useful in evaluating groundwater resources and solute transport in tectonically active areas. From this study, we can conclude that the tidal analysis of the groundwater levels can be used for reliable and accurate estimation of the aquifer’s properties and for determining their temporal variation.