1. Introduction

The frameworks of modern day characterization approaches are probabilistic and notoriously unreliable in application owing to inherent heterogeneity in fractured aquifers, scaling and uncertainty associated with their hydraulic parameter values; however, the cost dimensions of their hydraulic data requirements are not commensurate. Theories such as fractal transport are founded on knowledge of the probabilistic properties of the flow media [1] and [2]. As a potential measure of circumventing the high cost, [3] and [4] had conjectured that cheaper, minimally invasive and technically robust geophysical data ought to supplement hydraulic data. References [1] and [5] proposed a means of modeling groundwater flow dynamics based on crossover statistics of groundwater head fluctuations. It was observed that the multi-fractionality exhibited, was conditioned by other factors, scale-diversity in material-level heterogeneity, whose descriptions require robust and systematic quantification of multi-fractal behavior of the hydraulic conductivity, K field. Advances made recently on the application of metrics obtained via radar (GPR) technology [6] [7], Nuclear magnetic resonance, NMR [8] [9] [10] and seismic methods [11] have shown to correlate reasonably well with hydraulic conductivity. Despite the inroads made by various techniques, the convolution of poorly resolved computations of the parameter K persists.

In this paper, we advance theoretically and validate empirically a more deterministic electrical–hydraulic (eh) conductivity model in a fractured aquifer. The theoretical footing of the eh correlation rests on the stochastic models by [12] and [13] (hereinafter conveniently referred to as the model of Bernabé et al. and the BLM model respectively). The preferred “universal” power laws incorporated scale homogenization and uncertainty quantification; aspects which are conspicuously lacking in most probabilistic approaches. The BLM is a permeability k, model that quantifies the effects of pore network connectivity via the average coordination number z (i.e. the mean numbers of throats and pipes linked to a nodal pore), the hydraulic radius, r_H and many others. The incorporation of z in k, as a substitute to more expensive direct numerical simulation commonly applicable to macroscopic properties in the fine-scale [14], further restricts effects of model parametric uncertainty onto model outputs. The methodology is thus viable especially for clay-contaminated fractured aquifers.

The relation of Bernabé et al. attempts to identify micro-structural parameters influencing rock resistivity. It’s a permeability k, a formation factor F model; although it fails to unequivocally incorporate porosity $\varnothing$, it is congruent with Archie’s law in that the parameter z may be replaced by cementation factor exponent $m_{\varnothing }$ in Fa. Since determinations of hydraulic property by electric measurement is anchored essentially on the classical Archie’s relation, it seems practical that a micro-structural parameter must either increase or lose the connectedness of eh conduction space, for which both electrical and hydraulic properties of aquifer-level materials are contingent. By principle of conservation of connectedness, we formulate a hypothesis that whichever the structural eh conduction effect present in a 2-D sample, there is a dominant eh conduction phase. We then objectively validate the derived comprehensive eh conductivity relation in order to characterize a fractured rock aquifer wherein groundwater productivity is unsustainable.

2. Theory of Electrical—Hydraulic Conduction in Fractured Rocks

Minerals in rocks conduct electrical current through interstitial fluids present in the rock pore–field scale features. Reference [14] observed elevated insitu conductivities within a wide spectrum of 10⁻¹ to 10⁻⁸ mohm/m weighted against a narrow array of 10⁻¹² to 10⁻¹⁷ mohm/m in diverse rock materials. As a consequent, geo-electrical conductivity, c may be thought to be associated with fluid conductivity, C_f and pore-surface conductivity, C_s [15] [16]. The extended conceptualized models by [17] which comprises conductivity conditioned by: grain-matrix (C_gm), ion-exchange (C_ex), and the Maxwellian-effects (C_max) may be used to clear any doubts.

To begin with, the BLM model considered normal lattices bearing cylindrical pipes and circular cross-sections immersed in an impermeable matrix. To widen its range of applicability, they conceptualized 2-D and 3-D lattices. In the case of 2-D set-ups, a thickness was designated to the systems to compute k and $\varnothing$. By applying random (e.g. by changing pipe radii) network simulations and invoking the Kirchoff’s law under periodic boundary conditions, BLM concluded that k is a power function of z (except near the percolation threshold z_c) and pore radii distribution ${\sigma }_r$. Log-uniform distributions were applied since micro-structural parameters such as pore radii in rocks are robustly varied [18] [19]. Numerical solution sought at varying values of z and ${\sigma }_r$, yielded clustered average values of k and $\varnothing$. Their simplified relation was:

$k=\omega \frac{\pi }{8}{\left[\frac{r_H}{l}\right]}^2{\left(z-1.5\right)}^{\beta }{r}_H^2,$ (1)

where, $r_H$ is the hydraulic radius β and ω, are controls to the log-transformed measures of pore-size heterogeneity; and $\left[\frac{r_H}{l}\right]$, is the proportion of average effective pore radius to average pore parting transform (sometimes interpreted as longitudinal pore aspect ratio). The factor $\frac{\pi }{8}{\left[\frac{r_H}{l}\right]}^2$ is a consequence of Poiseuille formula.

As a way of enhancing the applicability of the model, the model of Bernabé et al. introduced the parameter, F to account for tortuosity in the regular networks of cylindrical conduits with circular sections immersed in an electrically insulated matrix. Unlike z in the BLM model, F is measurable through the Archie’s law. In order to simulate electrical conductivity through a network of pipes, a solution was sought by application of Kirchoff laws and the Poiseuille formula [12]. As for k in the BLM model, the model of Bernabé et al. studied the associations

of diverse permutations of $\varnothing$ and $\frac{1}{F}$ to $\left(z-z_c\right)$ and arrived at the same conclusion that $\frac{1}{F}$ obeyed a power law relationship stated as:

$\frac{1}{F}=w_F\pi {\left[\frac{r_H}{l}\right]}^2{\left(z-z_c\right)}^{\tau }.$ (2)

By combining Equations (1) and (2), the model of Bernabé et al. was simplified as:

$k=C{\left[\frac{r_H}{l}\right]}^{2\left(1-\propto \right)}{\left(\frac{1}{F}\right)}^{\propto }{r}_H^2,$ (3)

where C is a prefactor for circular pipes, $\propto ,\sigma$ are functions of scale variations in the model of Bernabé et al. This result may appear based on analytical relationship that electrical conduction and hydraulic flow replicas conform to equivalently recognized array of linear equations.

In spite of the advances to improve understanding of pore-level heterogeneity, it has been acknowledged that there does not exist a scale contingent universally accepted petrophysical models relating geophysical measurements to hydraulic parameters. In this study we perform modifications to scale-invariant stochastic relations above prior to applying them in a deterministic framework by using multi-scale data obtained from the Olbanita aquifer in Kenya Rift, for which so far no hydraulic characterization studies are available and sustainability of groundwater yield is uncertain.

3. Materials and Methods

3.1. Study Area

Olbanita aquifer occurs within the Lower Baringo Basin, Central Kenya Rift (Figure 1). The study area comprises a volcano-sedimentary sub-basin with an area of ≈6.25 sq. km roughly bound by latitudes 0˚07'S and 0˚10'N and longitudes 36˚05'E and 36˚15'S. The area has low relief (1650 - 1750 m asl.) and ASAL climate.

The aquifer is bisected at about four (4) places by roughly N-S normal faults a proxy probably showing portioning of the aquifer into hydraulic blocks. Within the blocks, flow could be factored by the rock matrix itself or a network of non-extensive fractures; herein conveniently referred to as Primary Fracture Network (PFN). Locally, the drainage could be from Menengai Caldera northwards; based on local piezo-metric levels. The N-S, NE-SW, and NW-SE form a Secondary Fracture Network (SFN) that provides channels diverting stream water underground at several places. Evidences based on borehole geo-logs indicate existence of a confined aquifer sandwiched between Wasagess flows (phonolites and trachytes) of the Rumuruti group and the Samburu basalts. Weathered and/or

fractured zones (comprised of trachytes tuffs and clays) in the volcanic rocks and sediments inter-bedded between volcanic rocks may also provide water in some localities.

3.2. Geo-Electric Resistivity Data

Vertical Electrical Resistivity (VES) datasets were the smallest scale of investigation adopted to characterize micro-structural parameters of the study area. The VES survey adopted the Schlumberger array configuration [20] using 19 profiles of 296 locations. Data inversion was conducted by using IP2Win software developed by Borland Int. and distributed by Geoscan M. Ltd, Moscow, Russia. The interactive-semi interactive concept of the computer code and inversion procedures used provided by IPI-2win. The VES data inversion output for each VES location was a graphical model defined by the layer resistivity ρ, and thickness h in plumb.

3.3. Aquifer Hydraulic Data

The aquifer hydraulic character in the Olbanita area was distinguished by confined weathered matrix interconnected and by non-extensive fractures (implying leaky conditions). Single well constant rate pumping test data collected recently by Rift Valley Water Services Board comprised the largest scale of investigation for this study. The test analysis has been, therefore, conducted by depending on Cooper-Jacob solution [21], stated as:

$Z_r=\frac{2.3Q}{4\pi T}{\log }_{10}\frac{2.25T{t}_0}{r_w^{2}S},$ (4)

where Z_r is the drawdown, Q is the pumping rate, T is the transmissivity, S is the storatility, t₀ is time interval since the start of pumping, and r_w is the bore radius (m).

3.4. Modifications to Pre-Existing Stochastic Models

Since detailed pore-level attributes are inconsequential individually but spatially varied bulk averages do, aquifer-level hydraulic determinations using stochasticity require up-scaling. Particle tracking studies by [22] suggest that fault conduits significantly influence results. K measurements by [23] imply a direct relation with density of conductive fracture zones and defects. However, Fracture density, hence matrix connectivity in volcanic rocks can increase with diminutive change in porosity [24], whereas at a local-scale higher porosity may slow-down the flow [25]. The Conduit Flow Process Model (CFPM) developed by [26] linked turbulent flow to a pipe-network observed in outcropping normal faults, implying viscous laminarity for freshly fractured rock aquifers [22]. Turbulent flow was also observed by [22] in the pipe network representing un-faulted rock matrix blocks away from fault zones in the equivalent porous medium characterized by Darcian fractures. Borehole geo-logs from the study area depict thick tuffaceous-clay layers. Constant rate aquifer tests also indicate low K values (in the order of magnitude of 10⁻¹) with a fluid exchange phenomenon in which block matrix flow aided by the PFN was identified. Further, the aggregate effect of water-rock interaction characterized by swelling of hydrophilic clays and consequent fracture space blockage implies that fractures undergo self-healing [27], gradually reducing k-values. Therefore, we conceptualized as most appropriate a turbulent flow regime (within both the partially closed rough fractures and block matrix) for the entire flow domain investigated.

The influence of such rough partially clay-filled apertures would be described better by using non-steady-state solution to the Darcian equation of flow

represented by; $Q=\frac{w}{12}\frac{e^3}{\mu }\frac{\Delta P}{L}$, where w is the width of each bond as detailed in

the Virtual Medial Axes (VMA) for fractured rock of [28], e is taken as the mean aperture thickness for an ideal fracture, μ is the groundwater dynamic viscosity, L is the path length where flow occurs and P is the pressure drop along the path. The hydraulic conductance K of ideal fractures would be distributed depending

on the cubic law stated as; $K=\frac{w}{12}\frac{e^3}{\mu }$; in which the more fundamental and intrinsic parameter k, which defines the pore geometry within the flow domain is represented by $k=\frac{e^2}{12}$. For the field hydraulic conductivity, K in Equation (3) can be re-written as:

$K=C\frac{we}{\mu }{\left[\frac{r_H}{l}\right]}^{2\left(1-\propto \right)}{\left(\frac{1}{F}\right)}^{\propto }{r}_H^2.$ (5)

Since geologic matrix hydraulic conductivity is based on aperture/pore conductance, the transmissivity, T for any aquifer interval is computed by using a geometric mean of aperture thicknesses a_c as distributed within the limits of the specified network interval. Therefore, we consider the dependence of field-scale effects on hydraulic apertures with a mean diameter, a_c, in the fracture-matrix network. The

network transmissivity would then be expressed as; $T=a_{c}K=a_{c}k\frac{we}{\mu }$. Unlike e,

a_c is a field-scale property which accounts for hydraulic channeling and the net fracture surface area [29]. Additionally, conventional engineering hydrogeology practice [12] [13] and [30] recommends the use, in place of a_c, a measurable

parameter such as hydraulic radius r_H (e.g. $r_H=\frac{2V_p}{A_p}$, where V_p and A_p are solid volume and exposed solid-fluid interfacial area of the fracture-pore matrix, respectively. Hydraulic radius is occasionally defined as; $\frac{V_p}{A_p}=\frac{r_H}{2}$ as a scale-invariant effect over a suite of aperture-networks with divergent pore radii distributions.

Rock formation at preliminary stages such as weathering and transport abrasion reorganize microstructural attributes proportionally [31] but randomly hence the expected variation in pore radius distributions and groundwater saturation. Locally, consequent distributions in electrical conductance depict broad and even skewed configurations [32]. Since the occurrences of micro-structural heterogeneities in rocks are strongly skewed, we consider the adequacy of log-uniform distributions of apertures, fractures, and so forth. The relation in Equation (5) becomes:

$\ln \left(T\right)=\ln \left[C\frac{we}{\mu }{\left[\frac{r_H}{l}\right]}^{2\left(1-\propto \right)}\right]-\frac{\propto }{r_H^{-3}}\ln F$ (6)

Equation (6) is, therefore, the modified aquifer-scale model which relates a borehole-scale parameter T, and a micro-scale parameter F. The equation is a linear function of the form:

$Y=a+bX,$ (7)

in which the coefficients a and b depicting the intercept and gradient between T and F, respectively, are stated as:

$a=\ln \left[C\frac{we}{\mu }{\left[\frac{r_H}{l}\right]}^{2\left(1-\propto \right)}\right]$ (8)

$b=-\frac{\propto }{r_H^{-3}}$ (9)

To achieve dimensional logic, the two coefficients are in units of m²/day, while F, being a ratio of resistivity values, is dimensionless. Since Equation (6) depicts field-scale parameters influenced by connectivity and heterogeneity of the flow medium, at each aquifer test site, natural logarithm of T was plotted versus natural logarithm of F derived from corresponding VES data.

Because of the flexibility of the theoretical eh network model (that is, it can accommodate any network heterogeneity), a bond shrinkage model (such as [31] model) was incorporated into the modified equation [15] [16] [17]. Accordingly,

the relationships between $\rho$ and k with $\varnothing$ fraction will be such that $\frac{1}{\rho }$ is expressed by Archie’s law; stated as $\frac{1}{\rho }\alpha \frac{1}{{\rho }_f}{\varnothing }^{\left(m_{\varnothing }\right)}$ whereas k is expressed by a Kozen-type power relation; stated as $k\alpha \frac{{\varnothing }^{\left(m_k\right)}}{S^2}$ wherein the power $m_{\varnothing }$ is

Archie’s exponent also known as cementation factor and S represents the specific surface area. Comparing these two relations, at macroscopic scale, $\varnothing$ is common in the numerator with positive exponents implying a positive eh correlation. Considerations of these proportionalities, e.g. [15-17] showed that:

$\ln T=\ln \left(\frac{B{a}_c\delta g}{\mu A}{S}^{-2}\right)-\frac{m_k}{m_{\varnothing }}\ln F_a^{-1},$ (10)

where in A and B are the pertinent proportionality constants $m_k$ and $m_{\varnothing }$ showing the different functions of bond shrinkage process; expressed accordingly

as; $m_k=\frac{2\ln x^2}{x^2-1}+\frac{2}{x+1}>0$ and $m_{\varnothing }=\frac{\ln x^2}{x^2-1}>0$; $0<x<1$, where, empirical

studies frequently report $m_{\varnothing }$ in the range of 1 - 2, and x is the bond shrinkage operator by which the pore radii of throats and conduits are adjusted [15]. The factors; $m_{\varnothing },m_k$ (hence the slope) must be positive since the pore-volume eh conduction envisaged by the bond shrinkage model depicts exactly this. To achieve structural synchrony between the modified models (Equations (6) and (10)) with

the empirical model of the Olbanita aquifer, its slope, $\frac{\propto }{r_H^{-3}}=0.537=\frac{m_k}{m_{\varnothing }}$ was expressed in relation to x as:

$b=\left[\frac{2\ln x^2}{x^2-1}+\frac{2}{x+1}\right]\left[\frac{\ln x^2}{x^2-1}\right]=\mathrm{0.537.}$ (11)

which further reduced to:

$2x-2=\ln x.$ (12)

Equation (12) was split into two relations stated as:

$y=2x-2,$ (13)

and,

$y=3\ln x.$ (14)

A numerical solution of Equations 13 and 14 wherein f(x) = y, y = 2x − 1 and y = 3 lnx for successive values of 0 < x > 1 was sought. When plotted simultaneously on one chart scale, the intersection of the two graphs yields the solution to the equations (Figure 2). Accordingly, the two functions intersect at x = 0.9 for values of $0<x<1$. For that reason, the cementation exponent $m_{\varnothing }$, defined

by $m_{\varnothing }=\frac{\ln x^2}{x^2-1}>0$; $0<x<1$, was 1.1.

The cementation exponent obtained was used in the computation of porosity values at each VES station by correct substitution in the proportionality relations given that EC of groundwater was observed at all the borehole sites. Field-wide hydraulic conductivity values were also obtained by using; $T=a_{c}K$, where fracture aperture $a_c$ is substituted by the more geologically plausible field-scale aquifer thickness $b_e\approx r_H$, obtained by taking the depth value of low-medium resistivity layer of the VES output model constrained by using borehole geo-log data. By re-stipulating that field-level eh relations are depended on mean sample $\varnothing$, we substituted the bulk resistivity alongside pore water resistivity into the

theoretically established Archie’s relation $\frac{\rho }{{\rho }_f}\approx \frac{1}{F}\propto {\varnothing }^{\left(m_{\varnothing }\right)}$ to obtain the resistive

component, F and porosity fraction, $\varnothing$ at each VES location. The macroscopic variation in K and F were then illustrated since the bond shrinkage model shows that as steps in the factor x get large, logarithmic ratios of the most probable mean ${\varnothing }_n$ evolve into the Kozen-type power equation and Archie’s law. (e.g. in [15]; after the work of [31].

4. Results and Discussion

4.1. Hydraulic Parameters

Aquifer storativity, flow and borehole yields in the area are found to be structurally controlled by the pervasive but variable array of fault-fracture traces. Groundwater BH 2 BH 5, BH 6, and BH 7 are intersected by major faults (Figure 1) and are the most productive in the area (Table 1). The Storativity, S-hydraulic conductivity, K relation depicts a direct correlation (highest in BH 6 and lowest in BH 1). BH 1 (the least productive) and BH 8 (dried-up borehole approx. 400 meters north of BH 7) are not intersected by any fault structure. The mean values of discharge recorded were highest at BH 6 (126.5 m³/hr) in the north western periphery of the study, whereas the lowest discharge boreholes, BH 7 and BH 7A with mean discharges of 57.5 m³/hr and 14.1 m³/hr, respectively, occurs in the eastern part near water plant. Groundwater boreholes located in the central portion (for example BH 1, BH 3 and BH 4) recorded average specific yields (that is, in the range of 1.48 - 1.86 m²/hr). BH 2, though located in this zone is one of the most productive probably due to deep weathering associated with a fault zone as depicted from its geo-log.

However, distinction of effective porosity values, hence permeability differentia in each continua becomes a crucial factor [22]. The intra-aquifer system fluid interaction observed in the data set was mainly block matrix flow to the pumping well via non-extensive fractures with the exception of BH 6. At this scale, pumping test data indicate that flow towards BH 6 is mainly controlled by SFN which is probably comprised of two sets: 1) the horizontal set responsible for sub-horizontal k, and 2) the less permeable sub-vertical set occurring within adjacent blocks that insures matrix connectivity decreasing drawdown at pumping.

Conversely, elsewhere the PFN influences the adjacent block matrix only at a decimeter scale, thereby increasing k and S of the block (Table 1).

4.2. The Calibrated Aquifer Model

The algorithm in Equation (6) was satisfactorily calibrated by comparison with more experimental data from the study area (Table 2). By applying Archie’s law, which has its own theoretical pedigree to data, it can be construed that the index F is less than unity, a clear indication that the rock matrix is less resistive than pore fluid. The resultant graph gives a positive gradient (Figure 3) indicating that groundwater flow is through non-extensive fracture-pore surfaces as opposed to flow through pore-volumes.

The observed linear regression correlation was:

$\ln T=0.537\left(\ln Fa\right)+3.695$ (15)

Equation (15) is the subsequent calibrated eh conductivity model approximated for Olbanita aquifer (Figure 3). The model is less rugged with reference

to a coefficient of correlation at 75.0% (percent) but we rest contented with the magnitude of precision since the transmissivity—resistivity factor correlation is “fit for purpose”.

The linearity of this graph was excellent with detailed resistivity data demonstrating power law relationships between field-based hydraulic parameters and micro-geometry based electrical parameters. By using numerical procedures described earlier, the study computed the cementation factor, m for the study area as $m_{\varnothing }=1.1$. However, the factor $\alpha$ in equation 19 must, therefore, remains negative in consistent with the positive slope observed in Equation (15), which depicts a negative eh conduction relation contrary to the bond shrinkage model envisaged by [31]. The low cementation exponent is generally indicative of clayey-contaminated fracture surface conduction characterized by cation exchange phenomena and/or fresh groundwater environment [33]. Rock fracturing, preferred flow paths along fracture-fissures and subsequent weathering may enhance contamination at these sites by alumino-silicates (clays and zeolites). Rapid rates of mineral dissolution are associated with areas where defects intersected the topography and vice-versa [34], causing modifications to pore throat aspects such as pore radii and pipe length (and hence porosity and permeability) variedly at different field sites [18]. Associated clay mineralization reduced the parameter z (hence low connectivity of pore volume networks) which obscures pore-volume electrolytic flow paths for continuous phase conductance, but favors fluid-solid interface (surface) conductance. The calibrated numerical model with deterministic range of borehole-scale conductivities and micro-scale formation factors produced an approximation of equivalent parameters at a higher scale. Geo-electric models of $\varnothing$, T, and K at each VES location are shown (Table 3). Apparently, at the scale of well-test (Table 1), K values observed are at least several orders of magnitude higher than the model values (Table 3), reflecting scaling effects. As [23] suggest, influences on small-scale measurements such VES, could probably emerge depended on the occurrence of sub-vertical artifacts (heterogeneity) and bias.

In this work, Archie’s law, and thus the positive eh correlation somehow fails to hold; probably signifying a pore surface dominant eh conductivity or freshwater saturated aquifer environment. An inverse eh conduction is probable in a pore surface Ap dominated conduction as opposed to a pore volume Vp environment [15]. The integrated effects of argillaceous materials (most commonly clayey minerals) coating pore surfaces lenders pore fluid conductivity too insignificant compared to pore surface conduction [35]. Presence of clay mineralogy may be attributed to alteration of feldspar minerals in volcanic tuffs clearly shown on most, if not all borehole geologic logs of the area. Clays by their inherent lattice structure provide negatively charged surfaces. Also, pH-conditioned adsorption reactions induce negative charges on clay surfaces. Whereas the former is widely considered permanent negative charge inscription, the latter (though less important) involves pore water dissociation whereby the hydroxyl

(OH⁻) are absorbed by clays and the H⁺ increases acidity of pore waters.

The consequent overall negative pore surface charge electrochemically bind positive charged ions initially dissolved in pore waters to pore walls, thereby forming an electrical double layer, which in turn acts as a surface for electro-migration of charges. By considering conductivity via the electrical double

layer as being dominant, and with the application of $F=\frac{c_f}{c}=\alpha \frac{1}{{\varnothing }^{\left(-m_{\varnothing }\right)}}$ and $k\alpha \frac{{\varnothing }^{\left(m_k\right)}}{S^2}$, specific surface area occurs in the numerator of the former (expressed

as pore network fluid conductivity, $c_f$ ) and $S^2$ in the denominator of the latter. The observation approves the inverse (1/F, K) correlation depicted in Figure 4.

A negative linear eh conductivity relation, therefore, arise which consequently implies a similar functional relation on a bi-log scale for this aquifer environment. Although it sounds awkward by using only a few data points in the high normalized electrical conductivity region, we remain unaware of any cases of equal dominance of two non-zero conductance phases in eh conduction space, even in multiphase mixtures. Nevertheless, the apparent realization of this model holds our hope that further works (with additional drill-sites) in the aquifer might 1) hold exactly or approximately true for Equation (15) or further approve the inverse (1/F, K) relation or 2) provide a relation that will further improve the model.

Meanwhile, at pore network scale, the Olbanita would be a case typical of low-salinity and/or clay-dominated aquifer that is characterized mainly by cation-exchange (CEC) phenomenon. Studies by [35] and [36] demonstrated that c_fαCEC. Under such scenarios, [15] urged that Archie-type and hence the positive eh conduction correlation is persistently restrained, setting forth a “less familiar” “Archie’s freshwater-clay counterpart”.

Results of this study agree reasonably well with hydraulic studies conducted by [37] [38] [39] [40] to investigate relations between pore’s internal surface area for each unit mass (m²/g) and CEC for freshwater (low EC) aquifers. Their laboratory samples were graded by grain size with a consequential finer grade in each successive sort. It was found that with each grain size reduction, electrical conductivity increased whereas K decreased. We also noticed that a transition from pore-volume to pore-surface conduction occurred after grain size decreased. In consistent with published literature, therefore, we concur that the finer the aquifer material, the greater the pore-surface area for electrical conduction through cation exchange (as in Equation (2)), and the higher the viscous drag for hydraulic flow (and consequent low K values) as depicted by Equation (6).

4.3. Conclusions

Adjustments have been made to the micro-level theoretical relations Bernabé et al. (2010, 2011) enhanced their applicability to field scale scenarios. From the widely known cubic law that controls flow through rock fractures, we expanded existing stochastic models by accounting for fracture wall roughness as well as inertia effects due to turbulent fluid in both the partially clay-filled fractures and block matrix. Specifically, the study modeled theoretically an association between field T and pore-level F. The network-level fracture-matrix heterogeneity and pore-level connectivity were found to be related by a negative linear function,

stated as; $\ln \left(T\right)=\ln \left[C\frac{we}{\mu }{\left[\frac{r_H}{l}\right]}^{2\left(1-\propto \right)}\right]-\frac{\propto }{r_H^3}\ln F$, in which the intercept is chiefly

a contingent of intrinsic properties of the eh conduction space and fluid property, whereas the gradient is based on the bond shrinkage factor x (for 0 1) of Wong et al. (1984) whereas x is determined by the mean eh conductive hydraulic radii over fracture-pore surface interface. The negative slope implies a theoretical inverse relationship between T and F as postulated in the pore-volume eh conduction case of Wong et al. (1984). By calibrating the modified model with transmissivity and formation resistivity data of the Olbanita aquifer, a positive linear relation of the theory, stated as $\ln T=0.537\left(\ln Fa\right)+3.695$, at R ² = 75.1% was obtained. By applying numerical procedures to the shrinkage model, the cementation factor, $m_{\varnothing }$ value for the aquifer material was 1.1. The empirical relation along with the low cementation exponent implies that T increases along with F, a case typical of eh conduction through pore-surface-networks contrary to pore-volume-networks of the shrinkage model case. By considering eh conductivity via the clay-surface electrical double layer as being dominant, we observed an inverse (1/F, K) correlation for this aquifer environment. The results of this study are in consistent with clay surface dominated eh conduction and/or fresh groundwater environment. Although proved robustly with further research, this is a valid hypothesis based on data deficiencies reported. Borehole geo-logs and aquifer-scale defects (fractures) provided further evidence of a weathering loci and consequent clay occurrences. The heterogeneous fractures provide preferred hydrologic paths, enhance rock dissolution and clay mineralization which consequently increase tortuosity and present significant complexity of electrolytic paths via pore-volumes. In principle, we conclude that fracture-controlled rock dissolution adjusts pore-level attributes and network connectivity; which in turn control k and $\varnothing$. We incontrovertibly highlight that the preferred eh conduction paths for the Olbanita aquifer were through interstitial solid-fluid interface (surface) cation exchange conductance. Though the modified relation is supported by many published fields and laboratory experiments for fractured-rock aquifers cited in this work, the calibrated relation may not be invariably applied to similar settings. Therefore, eh conduction modeling is a very reliable, less invasive, scalable and cost-effective method of hydraulic characterization in heterogeneous/fractured rock aquifers.

Acknowledgements

This study was funded by the Kenya Government through the National Research Fund http://researchfund.go.ke/. The original manuscript was also reviewed by an anonymous reviewer. The datasets used and generated in this study are available upon request.