Influence of Pore Size Distribution on the Electrokinetic Coupling Coefficient in Two-Phase Flow Conditions

: Streaming potential is a promising method for a variety of hydrogeophysical applications, including the characterisation of the critical zone, contaminant transport or saline intrusion. A simple bundle of capillary tubes model that accounts for realistic pore and pore throat size distribution of porous rocks is presented in this paper to simulate the electrokinetic coupling coefficient and compared with previously published models. In contrast to previous studies, the non-monotonic pore size distribution function used in our model relies on experimental data for Berea sandstone samples. In our approach, we combined this explicit capillary size distribution with the alternating radius of each capillary tube to mimic pores and pore throats of real rocks. The simulation results obtained with our model predicts water saturation dependence of the relative electrokinetic coupling coefficient more accurately compared with previous studies. Compared with previous studies, our simulation results demonstrate that the relative coupling coefficient remains stable at higher water saturations but vanishes to zero more rapidly as water saturation approaches the irreducible value. This prediction is consistent with the published experimental data. Moreover, our model was more accurate compared with previously published studies in computing the true irreducible water saturation relative to the value reported in an experimental study on a Berea sandstone sample saturated with tap water and liquid CO 2 . Further modifications, including explicit modelling of the capillary trapping of the non-wetting phase, are required to improve the accuracy of


Introduction
Groundwater provides the main source of water for human consumption worldwide, and it is also critically important for agriculture [1].Excessive abstraction of potable water from aquifers may cause exhaustion of the resource, especially during periods of reduced recharge.Moreover, poorly managed pumping schedules could potentially lead to contamination of boreholes, and in coastal aquifers, to the intrusion of saltwater [2].To improve aquifer management requires accurate characterisation and modelling of subsurface water flow using available methods.Recently, geophysical methods have been gaining a place of choice in the hydrogeologist toolbox, yielding the emergence of hydrogeophysics [3].Indeed, geophysical methods are non-intrusive and can be used for monitoring.However, such methods are not yet routinely deployed, and their implementation, as well as interpretation of measured signals, still require improvements.Thus, there is a real need for a robust hydrogeophysical tool to monitor water flow in hydrosystems.One of the main challenges is to characterise the effects of partial water saturation in the upper part of the critical zone, i.e., the vadose zone: the near-surface compartment where water and gas (normally air) are present.
The self-potential method (SPM) is a passive geophysical method that consists of measuring the naturally occurring electrical field (e.g., [4,5]).Among other applications, it has been shown to be a promising hydrogeophysical method to predict saltwater intrusion in coastal aquifers [6], to monitor rain or saline tracer infiltration in soils [7], to sense the root-water uptake linked to tree transpiration [8] or to characterise water flow through fractured systems [9].SPM is a non-intrusive, cheap and logistically simple method that can be used for long-term monitoring (e.g., [10]).It does not require large areas to lay out surface electrode arrays, as it can be deployed in boreholes.Moreover, the method is useful in relating the measured self-potential (SP) signal to aquifer heterogeneities, and the SP response also correlates with water saturation (Sw).
SP arises in subsurface settings saturated with water in response to pressure, concentration and temperature gradients, all of which may occur naturally in freshwater aquifers (e.g., [11]).Depending on the type of the thermodynamic driving force, SP can be categorised as electrokinetic or streaming (EK), electrochemical (EC) or thermoelectric (TE) potential if they are caused by pressure, concentration or temperature gradients, respectively (e.g., [4,12,13]).EK acts to maintain electroneutrality when an excess electrical charge at the rock-water interface is dragged by the flow.The excess charge that develops at the rock-water interface originates from the establishment of the electrical double layer (EDL) in response to the mineral surface charge that attracts ions of opposite polarity (counterions) from the bulk water (e.g., [14]).These counter-ions populate the diffuse part of EDL and are mobilised by the flow beyond the so-called slip plane, thus giving rise to EK potential (e.g., [15]).The electric potential at the slip plane is termed zeta potential, and it is a key petrophysical property that characterises rock-water interactions.The zeta potential can be estimated from the so-called EK coupling coefficient (CEK), which is defined as a ratio of voltage to the pressure difference, and it can be routinely measured in the subsurface and laboratory (e.g., [16,17]).
Vinogradov and Jackson [18] experimentally demonstrated that CEK depends on Sw.To demonstrate the correlation, Vinogradov and Jackson [18] used the relative coupling coefficient (Cr) that relates CEK at partial saturation to that at Sw = 1.Previous empirical, analytical and numerical studies have tried to relate Cr(Sw) to various aquifer and fluid properties (e.g., [15,[19][20][21][22]).However, there is still no consensus in the literature, thus suggesting that although Cr(Sw) is a crucial parameter for flow characterisation, it is still poorly understood.Moreover, it appears that Cr(Sw) strongly depends on the porous medium considered in the study [22].
Measuring Cr(Sw) is a challenging task, which cannot be carried out for all possible permutations of minerals and fluids.Moreover, the behaviour of Cr with water saturation is not always predictable or easy to interpret from experimental observations.Therefore, it is essential to develop a model that explains the expected Cr(Sw) and, more importantly, is capable of predicting the correlation between Cr and Sw.The current modelling techniques of Cr consider four main approaches: (1) the pore-network modelling (e.g., [23,24]), which represents the pore space topology as an ensemble of larger voids (pores) connected by narrower channels (pore throats); (2) direct pore-scale modelling in which the governing equations are solved for a precisely reconstructed pore-geometry (e.g., [25]); (3) the Representative Elementary Volume modelling that represents a porous medium as a number of blocks (elements) so that the fluid properties are assigned to each block, and the governing equations are solved using finite difference (volume) or finite elements approach (e.g., [22]); (4) the Bundle Of Capillary Tubes (BOCT) modelling that uses a bundle of parallel tortuous capillaries as a representation of the pore space (e.g., [26,27]).
All the above approaches have certain advantages and drawbacks.The pore-network and direct pore-scale models require precise reconstruction of the pore space from SEM and/or micro-CT images, both of which are technologically challenging and time-consuming.Moreover, the produced images by SEM and micro-CT are of nm-to μm-scale and therefore, the image-based models are of the same scale, which in most cases cannot be representative of a larger scale experimental condition, and the models need to be up-scaled, thus losing their precision.Despite the fact that these models are efficient in solving the governing equations, image processing is computationally expensive and timeconsuming.Moreover, unlike the direct pore-scale model, the pore-network models represent the pore space by spherical voids and circular (triangular, rectangular, etc.) ducts (e.g., [23]); hence, these models do not capture true pore space topology.Even though Jougnot et al. [24] showed that this true pore geometry is not necessarily needed to predict electrokinetic couplings in saturated conditions in the absence of surface electrical conductivity, the pore space topology plays a major role when considering the partially saturated conditions.
On the other hand, the REV approach can be applied for modelling flows on a much larger scale (up to kms), and its computational efficiency depends on the number of gridblocks used to model aquifers, so that more accurate representation of subsurface setting would require millions of grid-blocks and the simulation would become extremely timeconsuming.Moreover, such REV models take average fluid and aquifer properties across the size of a grid block, thus disregarding accurate physics of multi-phase flows.Simulated multi-phase flow using REV approach also relies on inaccurate, often empirically assumed, relative permeability and capillary pressure data (e.g., [22]).
The BOCT approach in modelling multi-phase flow in porous media is simple and captures and explicitly describes most of the petrophysical properties, such as water content (e.g., [28]), porosity and permeability (e.g., [29,30]), electrical conductivity (e.g., [31,32]) and thermal conductivity (e.g., [33]).However, the representation of the pore space is quite poor since pores are represented by capillary tubes that do not intersect and fluids flow in one direction.Such an approach is, of course, useful in modelling real mmto cm-scale coreflooding experiments, in which flows are 1-D in nature.Thus, the BOCT model is capable of providing useful insights into single-and multi-phase flows in porous media on a larger scale than pore-network or direct pore-scale models while still accurately capturing basic petrophysical properties from the fundamental hydrodynamic principles.There are, however, some limitations of the BOCT approach that include unrealistic capillary (pore) size distribution and constant capillary radius even for 1-D simulations, and which need to be addressed and improved.
One can utilise the approach of Jougnot et al. [27] to determine the EK coupling coefficient based on inferred capillary size distribution.In this study, the authors use two hydrodynamic functions (the so-called water retention and relative permeability function) to infer the size distribution of equivalent straight capillaries.In their study, the EK coupling coefficients predicted from their model are highly soil-or rock-specific.
Therefore, the aim of this study is to improve on published BOCT models by introducing a more realistic non-monotonic capillary size distribution obtained from direct pore size distribution measurements and alternating capillary radii.We report the effect of the new capillary size distribution on simulated petrophysical properties in comparison with previously published studies.Moreover, we report a case study and analysis of the application of the novel BOCT model to experimental data on partially saturated porous media reported in the case study.

Basic Definitions and Capillary Size Distribution
The classical BOCT model considers N capillaries, with each capillary permitted to have a different radius rc.The capillaries are confined within a model of length L and area A, as displayed in Figure 1.Each capillary within the model is defined by a length Lc, tortuosity = / and a cross-sectional area to flow = [26,34].In the alternating capillary radius model, the pore throat is modelled by a reduced radius , while the pore throat length is defined by , where is the pattern length, which is described in more detail below [35].
The bundle of capillary tubes model deviates from the true pore structure of porous geologic media, as it is formulated on the assumption that there are no points of intersection between capillaries and that all capillaries run parallel with the same orientation [26].It is widely acknowledged that the pore structure of geologic media is extremely complex, resulting in the effective porosity and absolute permeability of the media being highly dependent on the ratio of interconnected pore space [36].However, for the simplicity of modelling, the assumption of zero intersections ensures that the transport of mass and electrical charge through the model takes place in only one direction [26].Furthermore, in order to define an average tortuosity t for the model, it is assumed that the capillary radius rc and length Lc are independent of each other.
The previously published works of Jackson [21,26] and Soldi et al. [37] define the distribution of capillary radii throughout the model differently.Jackson [21,26] states the number of capillaries with radius between and + as ( ) such that the total number of capillary tubes is where N is the total number of capillaries, and and denote the minimum and maximum capillary radius within the model, respectively.The distribution of capillaries throughout the model is associated with the capillary radius by a function of the form [25].
where DJ is the normalization factor (constant, unitless) and 0 < m < ∞ (constant, unitless).The distribution function monotonically decreases at all values of m with the exception of m = 0, where the function produces a uniform distribution of capillary radii between and (Figure 3 in [26]).As the value of m increases, the frequency distribution is skewed more towards smaller radii values, similar to the skewed pore size distributions exhibited in many geologic porous media [26,38,39].
The frequency distribution of capillary radii presented by Soldi et al. [37] is somewhat different to that of Jackson [21,26].Soldi et al. [37] assume that the number of pores with radii greater than or equal to a specific capillary radius R follows a fractal law such that [28,29,35,40]: where 0 < ≤ ≤ < , DS represents the dimensionless pore fractal dimension (1 < < 2) and RREV is the radius of a cylindrical representative elementary volume (e.g., a cylindrical sample of porous geologic media) [37].If = , then = 1 and the cylindrical REV is entirely occupied by a single pore [35].Conversely, if = 0, there are an infinite number of pores within the REV [35].The cumulative pore size distribution (Equation ( 3)) is then differentiated with respect to to obtain an equivalent expression to Equation (2).The distribution of capillaries throughout the model is defined in terms of ( ), which represents the number of capillaries whose radii fall within the small range to + [37]: A significant difference between the distributions produced from Equations ( 2) and ( 4) stems from a criterion of the model proposed by Soldi et al. [37] that the largest capillary must always exist.Unlike the frequency distribution of Jackson [21,26] (Equation ( 2)), that of Soldi et al. [37] (Equation ( 4)), cannot be skewed and so the accuracy of choice for rmax greatly influences how representative the BOCT model is of the modelled porous geologic media.
Nonetheless, by the implementation of capillary tubes with varying radii, the work of Soldi et al. [35,37] enhanced the BOCT model.The addition of pore throats within the model, rather than the constant radius approach adopted by Jackson [21,26], provides a more realistic representation of the pore structure in geologic media, thus enabling explicit modelling of residual phase saturation and impact of wettability alteration.The inset of Figure 1 displays the pore geometry of a single capillary tube proposed by Soldi et al. [35].The radius at a particular point x along the length of a capillary is expressed as [35]: where a denotes the ratio in which the capillary radius is reduced, known as the radial factor (0 ≤ ≤ 1, dimensionless), hence represents the radius of the periodically distributed pore throats along the length of each capillary.It is assumed that each pore and its adjacent pore throat have a wavelength, or pattern length, , with the length of each capillary containing an integer number of wavelengths, .Subsequently, the length of a pore throat is expressed as , where represents the segment of which contains a pore throat, also known as the length factor ( 0 ≤ ≤ 1 , dimensionless), and = 0, 1, … -1 [35].
Although the bundle of capillary tubes model is not the most accurate representation of the pore space in geologic media, the published works of Jackson [21,26] and Soldi et al. [37] individually provide good insight into the applications and advantages of modelling porous media as a bundle of tortuous capillary tubes.However, neither Jackson [21,26] nor Soldi et al. [37] considered non-monotonic capillary size distribution.On the other hand, Jougnot et al. [27] considered non-monotonic capillary size distribution but used the constant capillary radius approach, and their capillary size distribution was indirectly inferred from other modelled properties rather than from actually measured one.Moreover, many published studies investigated such non-monotonic pore size distribution and its effects on capillary pressure [41], relative permeability [42] or hydraulic conductivity [43].However, to the best of the authors' knowledge, such capillary distribution functions have never been used in modelling coupled electro-hydrodynamic problems.In this study, we present a new BOCT model that considers a non-monotonic pore size distribution based on experimentally measured pore size distribution and alternating capillary radius.Moreover, our new model is constructed to match the reported porosity and permeability values of a real sandstone sample and cross-compare the modelled hydro-and electrodynamic properties between the three types of capillary size distribution (termed Jackson, Soldi, New) and both constant and alternating capillary radius approaches.It should be noted that although our BOCT model was developed using available information on pore size distribution for a given rock sample, unlike in pore-network or direct pore-scale models, such information is not essential for BOCT models as capillary (pore) size distribution can be obtained from matching sample's porosity and permeability.

Minimum and Maximum Capillary Radius
The Berea sandstone core sample, as discussed in Moore et al. [44], was modelled as a bundle of tortuous capillary tubes, assuming the system was water wet.Prior to implementing any capillary size distribution (CSD) function, a review of relevant literature surrounding Berea sandstone pore size distribution was conducted in order to determine the most appropriate minimum and maximum capillary radius to be used in the model.Table 1 summarises the six key research papers consulted during the decision-making process, and Figure 2a,b present the pore size distributions from Shi et al. [39] and Li and Horne [38], respectively.The values of pore radius and throat radius given in Table 1 are either explicitly stated within the respective research paper as a single value or have been taken as the peak value from a pore size distribution curve.The findings in the paper of Li and Horne [38] are presented in terms of pore throat diameter.Therefore, a rough assumption was made that the pore throat diameter was approximately equal to the pore radius, assuming a radial factor (ratio between throat and capillary radius) of 0.5, and so the pore radius has been reported as 5 μm in Table 1.Furthermore, Minagawa et al. [46] present a pore and throat size distribution as one figure, the peak of which sits at approximately 14 μm.Therefore, it has been assumed that this figure represents the average of all pore and throat radii, and so 14 μm has been presented in Table 1 as both the pore radius and the throat radius.The radius values from Ott et al. [47] and Thomson et al. [48] were determined to be outliers as they were, respectively, significantly larger and smaller than other values from the literature.
Initial values of rmin and rmax were considered by consulting the pore radius and throat data presented in Table 1 and Figure 2 and chosen to be 5 and 60 μm, respectively.The minimum value of 5 μm was selected due to the negligible pore volume of capillaries below this value, as depicted in Figure 2a, which will subsequently only be occupied by the wetting phase at the irreducible water saturation and will not contribute to flow within the core sample-therefore, having no impact on the permeability of the model.The minimum rmax value considered was 60 μm, again selected based on the maximum radius presented in Figure 2a.However, the simulation procedure was conducted using rmax of both 60 μm and 100 μm to investigate the influence of rmax on the distribution functions.

Matching Porosity and Permeability to Berea Sandstone
The porosity and permeability of the BOCT model were calculated using Jackson (2010) (Equations ( 4) and ( 5), respectively).Three CSD functions were then investigated, with their constants and exponents adjusted accordingly, in order to match the porosity of the model to the porosity of the Berea sandstone sample reported as 18.75% in Moore et al. [44].Since no permeability was provided by Moore et al. [44], the permeability of the Berea core sample, and therefore BOCT model, was estimated from the literature data obtained from experimental procedures (Table 1).This resulted in a permeability range of 10s to low 100s mD being considered acceptable for a sample with 18.75% porosity.
Accurately predicting the tortuosity of the Berea sample was key in matching all three CSD functions, as tortuosity is required to calculate both porosity and absolute permeability using both Jackson and New approaches (note the porosity and permeability calculation using New CSD relies on Jackson [21] Equations ( 4) and ( 5)).To compute the sample porosity using the Soldi CSD also requires tortuosity (Equation ( 5) in [37]).The tortuosity of the model was calculated using studies of Civan (2007, Equation ( 3)-( 14)) and Jougnot et al. ( [15], Equation 4.29), in which F represents the formation factor of the Berea sandstone, which was explicitly stated as 24 in Moore et al. [44].The Civan [49] and Jougnot et al. [15] approaches yielded tortuosity of 4.33 and 2.12, respectively.An additional calculation method outlined by Tiab and Donaldson [50] returned a tortuosity of 10.77, whilst Zecca et al. [51] recorded tortuosities of 4.3 and 4.5 for Berea sandstone using water and methane at 3 MPa, respectively.Attia et al. [52] measured tortuosities in the range of 3.78-4.71for Berea sandstone, with an average value of 4.16.Therefore, range of values was investigated to determine the most realistic tortuosity and value of 4.5 was selected for a sample with porosity of 18.75%.In this study, we assumed the tortuosity to be independent of phase saturation, consistent with published studies [53,54].

Constant Capillary Radius
The first CSD examined was that proposed by Jackson [26] using Equation (2) and constant capillary radius.When manipulating the m and DJ parameters of the Jackson CSD, it was imperative to decide on whether the total number of capillaries was kept constant or the total volume of capillaries was constant.As and were predetermined based on literature values, in order to accurately replicate the sample examined by Moore et al. [44] for which the absolute permeability was not reported, the values of m and DJ were optimised until the porosity of the BOCT model matched that of the Berea core sample (18.5%), therefore ensuring the total volume of capillaries remained constant.Furthermore, it was crucial to ensure that the CSD resembled the pore size distributions depicted in Figure 2a (the right-hand side at pore radii greater than 6 μm) whilst still ensuring a realistic permeability for the model consistent with literature values (Table 1).It was, therefore, important to first choose the exponent m that gave the general shape of the distribution required and then to adjust DJ to achieve a porosity of 18.75%.
The second capillary size distribution examined was Soldi CSD using Equation ( 4).As all values of had already been incrementally distributed between and and RREV was a fixed value of 12,500 μm taken from Moore et al. [44], the only parameter which could be adjusted to match the desired porosity was the fractal dimension, DS.Therefore, the fractal dimension was manipulated until the porosity of the BOCT model matched that of the Berea sandstone sample.
Following the review of the literature surrounding Berea pore size distribution, the findings of Shi et al. [39] were determined to be the most similar to the core plug examined by Moore et al. [44].This conclusion was drawn as the porosity reported in both papers was the same, and the permeability reported by Shi et al. [39] fell nicely within the range of permeability expected of Berea sandstone.Furthermore, Shi et al. [39] presented a detailed pore size distribution as displayed in Figure 2a, and for these reasons, the New CSD for the BOCT model was predominantly based on these findings.
The desired shape of the New CSD was replicated by modifying Jackson [26] distribution (Equation ( 2)) and invoking three intervals such that to capture the peak pore radius observed in the literature: The limits of each interval were determined by calculating the average pore size from the data in Table 1, which resulted in an average peak pore radius of 11 μm, thus informing on the location and the width of the peak.As before, and represent the normalisation factors of the first and third interval, respectively, with m1 and m2 being the respective skewing constants.
Table 2 displays the final values required for the parameters of all three distribution functions for the BOCT model with constant capillary radius, each resulting in a porosity of 18.75%.Bearing in mind the typically expected permeability of Berea sandstone from literature (Table 1), it was concluded that the values in Table 2 resulted in the most accurate and realistic pore size distributions and permeabilities for the model with 18.75% porosity, with Table 3 providing the resulting permeability and total number of capillaries within the model.It is clear that the permeability obtained with Soldi CSD using Equation ( 4) is outside the expected range of the literature.This is caused by the inability to skew the Soldi CSD function towards smaller capillaries due to the criterion of the model that the largest capillary must always exist. Figure 3 presents a comparison of each distribution function using Equations ( 2), ( 4) and ( 6) with a maximum capillary radius of 60 μm (Figure 3a) and 100 μm (Figure 3b).It can be seen from Figure 3 and Table 3 that the New CSD provides a more realistic size distribution as it is non-monotonic, as expected from real rocks, and it results in significantly fewer capillaries with radius less than 10 μm compared with Soldi and Jackson CSD.Both the Soldi and Jackson CSD require a large number of small capillaries to accommodate for the porosity of 18.75%, while these capillaries do not contribute significantly to flow and hence the permeability.Moreover, increasing rmax to 100 μm results in nearly 3-fold increase in the model permeability using the Soldi and Jackson CSD functions, while the permeability obtained with the New distribution remains fairly constant, which is consistent with real rock permeability behaviour if very few larger pores are added to (found in) the sample.

Alternating Capillary Radius
To model the CSD with alternating capillary radius requires determination of the radial (a) and length (c) factors in Equation ( 5).Examining micro-CT images and pore network modelling, such as the work of Sharqaway [55], allowed average a and c values to be estimated.Furthermore, the published work of Shehata et al. [56] determined the pore throat radius of Berea to be in the range of 2.5-25 μm, therefore permitting minimum and maximum radial factors of approximately 0.5 and 0.85, respectively.Upon sensitivity analysis of the a and c (not presented here) and consultation of Sharqaway [55] and Shehata et al. [56], a radial factor = 0.7 and length factor = 0.2 were decided to be the most representative of Berea sandstone whilst returning desired permeability values.

Fraction of capillaries, %
Fraction of capillaries, % Equations for porosity and permeability of the model were altered to account for the implementation of pore throats.The porosity of the model was calculated using an altered form of Equation (4) from Jackson [21] such that: where t is the previously defined tortuosity of the model, A is the cross-sectional area of the Berea sandstone sample and ( , ) is a dimensionless volume factor (same as Equation (6) in Soldi et al., [37]): ( , ) = ( − 1) + 1 Similarly, the permeability of the model was calculated using an altered form of Equation ( 5) from Jackson [21] such that: where ( , ) is a permeability factor that quantifies the reduction in volumetric flow rate due to the pore throats (same as Equation ( 18) in [37]): Implementing the radial and length factors into Equations ( 8) and ( 10) resulted in volume and permeability factors, fv and fk, of 0.898 and 0.612, respectively.
Table 4 displays the final values required for the parameters of all three distribution functions-Equations ( 2), ( 4) and ( 6)-for BOCT model with alternating capillary radius, each resulting in a porosity of 18.75%, with capillary radii distributed in increments ( ) of 0.1.The permeability and total number of capillaries within the bundle of capillary tubes model using each of the distribution functions is detailed in Table 5.The addition of pore throats resulted in a permeability reduction of approximately 83% with all distribution functions.Since the porosity of the model was reduced by the volume factor, the number of capillaries within each model also increased to achieve the desired porosity of 18.75% by 11.3% for both the Jackson and New distribution functions, and 12.9% and 13.7% for the Soldi distribution model with rmax of 60 μm and 100 μm, respectively.The shape of the capillary size distribution curves remained the same as before with the constant capillary radius model for all three distribution functions, i.e., the m, m1 and m2 exponents required for the Jackson and New distribution functions were untouched.The only aspect of the CSD that changed was the number of capillaries of each size (y-axis magnitude), as apparent in Figure 4.The change is due to the reduced volume of each capillary and thus the need for a greater number of capillaries within the model to achieve a porosity of 18.75%.

Multi-Phase Flow Simulation
The model was then extended to multi-phase flow by considering a wetting phase w (water) and an immiscible non-wetting phase nw (CO2 as in [44]).Based on experimental data [57], Berea sandstones exhibit strongly water-wet behaviour when saturated with liquid CO2 and high salinity brine.In the experiments conducted by Moore et al. [44] Berea sample was also saturated with liquid CO2, but tap water was used.Considering the thermodynamics of wettability [58], it is expected that water wetness should increase with decreasing salinity.Therefore, we assumed that our BOCT model is strongly water-wet, and we only consider this case.Moreover, an assumption was made that capillaries occupied by the non-wetting CO2 contain a thin immobile layer of water [21].This volumetrically insignificant layer of water was included in the model as it contributes to the surface electrical conductivity and regulates the development of an electrical double layer in the wetting phase (e.g., [59]).It was also assumed that the non-wetting CO2 is non-conductive, and therefore there is no electrical double layer in the non-wetting phase and associated with its surface electrical conductivity.
To determine water saturation, Sw, and the relative permeability to wetting phase, krw, of the model with constant capillary radius, we used the approach described by Jackson ( [21]; Equations ( 6) and ( 7), respectively).explicitly capturing the residual trapping of the non-wetting CO2, hence allowing the entire capillary to be occupied by either wetting or non-wetting phase.In reality, in some capillary tubes, depending on the capillary entry pressure, only pores (and usually not pore throats) should be occupied by the non-wetting phase via the so-called snap-off mechanism, which would have an impact on the relative permeability.However, it is clear from Figure 5 that the distribution of capillary radii throughout the model has a significant impact on the gradient of the krw curve at low irreducible water saturations, with each distribution producing notably different values of krw in the domain ≤ 0.8 when = 0.The model constructed using the New distribution function has larger relative permeabilities at lower water saturation, whilst the relative permeability of the Soldi model increases notably slower.As Swirr increases towards 0.4, the krw curves of each of the three distribution functions gradually move closer together and begin to overlap.These results suggest that the distribution of capillary radii throughout the BOCT model becomes less significant to the relative permeability of the model as the irreducible water saturation increases.d-f) of the BOCT model constructed using the Soldi (a,d), Jackson (b,e) and New (c,f) distribution functions with capillaries of alternating radius assuming a thick electrical double layer and significant surface conductivity with maximum capillary radius of 100 μm.

Relative Streaming Potential Coupling Coefficient
The relative streaming potential coupling coefficient Cr was investigated for all three distribution functions invoking both the thin and thick double-layer assumptions.It was found that alternating capillary radius did not have any noticeable effect on Cr for either CSD or whether the thin or thick double-layer assumption was invoked.This is explained by the inability of the current model to capture realistic residual trapping of the non-wetting phase, thus assuming that any capillary tube in its entirety is occupied either by water or by CO2 regardless of whether the corresponding capillary has a constant or alternating radius.Moreover, varying the capillary size distribution function had no noticeable effect on Cr when employing the limit of a thin double layer.This conclusion was drawn when considering both zero and non-zero ( is less than 10% of ) surface conductivities.Moore et al. [44] used tap water of 125 Ω • m electrical resistivity in their streaming potential measurements.Another experimental study on multi-phase streaming potential measurements [18] reported values of Cr as a function of Sw using liquid non-polar undecane and 0.01 mol/L NaCl solution of approximately 10 Ω • m resistivity.Consistent with previously reported results [60], surface conductivity becomes dominant when the bulk resistivity exceeds 10 Ω • m.Therefore, we apply the thick double-layer assumption to describe the water saturation dependence of Cr and compare the modelling results to the published experimental data.Figure 5d-f displays the relative coupling coefficient results for all three capillary size distribution functions investigated, assuming a thick electrical double layer, with model = 100 μm and capillaries of alternating radius at various Swirr values from 0 to 0.5.We do not present results obtained with constant capillary radii as they are essentially identical to those with alternating radii.The difference between the three capillary size distribution functions became apparent when analysing the Cr results considering a thick electrical double layer-particularly at small values of Swirr (between 0 and 0.2).The New distribution function presents the sharpest decrease in Cr as Swirr approaches 0 in comparison with other CSD functions, which mimics the results reported by Vinogradov and Jackson ([18]; Figure 3a,c) and simulated by Zhang et al. ( [22]; Figure 13a,c).The results presented in Figure 5 demonstrate that more accurate modelling of the capillary size distribution is capable of capturing the experimentally observed behaviour of the Cr(Sw).However, the Cr behaviour with Sw obtained using the Jackson CSD closely resembles the results of the New CSD.On the other hand, the Soldi CSD does not provide as steep of decrease in Cr at proximity to 0 of Sw, thus suggesting this approach is least suitable for simulating the electrokinetic properties of the type of rocks studied here due to the very large number of small capillaries in the fractal distribution.Moore et al. (2004) In order to simulate Cr of the Berea sandstone sample, the surface conductivity was determined to match the reported water electrical conductivity and that of the fully water-saturated rock sample of 3.03 × 10 −3 S/m [44].

Bundle of Capillary Tubes vs. Experimental Results of
To determine the surface conductivity of the Berea sample, MATLAB (Math-Works, Natick, MA, USA) scripts were written to calculate the effective conductivity and conductance of a single capillary, the total conductance of all capillaries within the model constructed using each of the three distribution functions, and finally, the conductivity of the water-saturated sample: where L and A are the length and area of the Berea sandstone sample, respectively.The results for each distribution function and rmax values are displayed in Table 6.As the surface conductivity is multiplied by the number of capillaries with a radius between and + , the distribution of capillary radii throughout the model is therefore a limiting factor of the surface conductivity and, as such, each distribution at each value of rmax resulted in different surface conductivities.Comparing the results of Table 6 with those of Table 5, it is suggested that the surface conductivity correlates with the permeability of the model in some way, as the surface conductivity increases with increased permeability.This makes sense as the decrease of permeability is due to an increasing number of small capillaries, hence an increase of specific surface area and therefore more surface with an electrical double layer generating this surface conductivity.Due to the formulation of the New distribution function and subsequent similarity in model permeability independent of the rmax value selected, the New distribution was the only distribution that resulted in similar surface conductivities for rmax of 60 μm and 100 μm.As the bulk conductivity of tap water within the pore space of the Berea sample was measured to be 8 × 10 −3 S/m (8 × 10 −9 S/μm) by Moore et al. [44], the bundle of capillary tubes model is dominated by the surface conductivity since it is three orders of magnitude larger than the bulk conductivity in capillaries of an order of 1 μm and one order of magnitude larger than the bulk conductivity in the largest tested capillary tubes of 100 μm (note that the bulk conductivity is multiplied by the term and the surface conductivity is multiplied by the term in Equation ( 16)).The high surface conductivities presented in Table 6 are indeed expected as the bulk electrolyte (tap water) salinity is at least 1 order of magnitude lower than the threshold salinity at which the surface conductivity becomes dominant in sandstones [60] and comparable to the corresponding value for sand packs [61].The streaming potential coupling coefficient of the fully water-saturated rock sample was calculated with respect to the zeta potential using: Invoking thin double layer and assuming the electric permittivity and viscosity of water to be 7.1 × 10 F/m and 8.9 × 10 Pa • s, respectively, with the aim of achieving the coupling coefficient of −300 mV/MPa as reported by Moore et al. [44].Although Jackson [21] and Linde [62] derived expressions to compute ( = 1) in the limit of a thick double layer, the proposed equations use the surface charge density, which is neither reported nor can be interpreted from Moore et al. [44], and for that reason, we used Equation (17).The bulk conductivity of the model was set to the measured value of the core sample (8 × 10 −3 S/m), and the zeta potential ζ was assumed to be −28 mV as calculated by Moore et al. [44].The coupling coefficient at = 1 was calculated for each distribution function with rmax values of 60 μm and 100 μm using the respective surface conductivities determined previously.
Unsurprisingly, the simulated coupling coefficient for the fully water-saturated rock sample under each of the distribution functions does not match the coupling coefficient of −300 mV/MPa reported by Moore et al. [44]-the simulated coupling coefficient is −62 mV/MPa, which is approximately five times smaller.However, in order to simulate liquid CO2 flooding of a tap-water-saturated sample and to match the relative coupling coefficient of the model to that of Moore et al. [44], it is not necessary to match ( = 1) of the model to that of Moore et al. [44].
The relative streaming potential coupling coefficient of the BOCT model was simulated under the thick double-layer assumption using Equation (15), with the aim of replicating the relative coupling coefficient of 0.1 reported by Moore et al. [44].The simulated relative coupling coefficient results are presented in Table 7 for models with = 100 μm, and which summarises the actual Swirr of the Berea sample for a range of possible values which could be reported by Moore et al. [44].Based on the Swirr values reported in the literature for Berea sandstone, Cr was initially simulated for Swirr in the range 0.2-0.5.Since according to our model, a non-zero coupling can only occur when water is flowing within the Berea sample, acknowledging that if water is still flowing, then the irreducible water saturation has not yet been reached.This is consistent with the conclusions drawn by Vinogradov and Jackson [18] and Zhang et al. [22].On the other hand, an experimental study published by Allègre et al. [63] reported a non-zero streaming potential coupling coefficient at irreducible water saturation and presented an empirical model to explain this behaviour.However, the model heavily relied on an abnormal relative streaming potential coupling coefficient of an order of 30 at partial water saturation.To the best of our knowledge, such large values of Cr have not been reported by any other research group, hence the empirical parameter introduced by Allègre et al. [63] that allowed Cr at partial saturation to be 200 greater than the corresponding value at Sw = 1 was considered to be irrelevant for our model describing experimental conditions when ( ) ≪ 1.The results in Table 7 are designed to be analysed such that if Moore et al. [44] report ( = 0.32) = 0.1, the true Swirr should be 0.21, 0.22 and 0.22 for the Soldi, Jackson and New CSD functions, respectively.Therefore, Table 7 effectively demonstrates the magnitude of the error within Swirr of Moore et al. [44], with the largest and the smallest discrepancies between the reported and actual Swirr being established with the Soldi and New distribution functions, respectively.However, at the reported Swirr of 0.3 and above, the actual Swirr of all distributions is approximately the same.
Following the comparison of the BOCT model results and experimental results of Moore et al. [44], this work confirms that the Moore et al. [44] paper must be corrected for the actual Swirr, since it is not possible to achieve a non-zero coupling coefficient at the irreducible water saturation due to the lack of a flowing electrolyte.Our model also suggests that a more accurate New CSD function should be used in BOCT models to capture the behaviour of Cr with water saturation.Moreover, it is suggested to explicitly implement the residual trapping of a non-wetting phase in the BOCT model with alternating capillary radius.

Conclusions
The bundle of capillary tubes model developed in this study to simulate the streaming potential coupling coefficient in partially saturated porous media is based on a realistic pore and pore throat size distribution of real rock samples.Compared with previously published studies, our model produces results that are more consistent with existing experimental data.


Unlike the previous bundle of parallel capillaries models, our approach to defining the pore and pore throat radii distribution is based on direct measurements, thus providing a more realistic description of porous rocks, in which pore and pore throat size distribution is non-monotonic;  Our model was tested using constant and alternating capillary radii, with the latter being invoked in order to distinguish between pores and pore throats.Despite the alternating capillary radii model's capability, we did not to explicitly model residual trapping of the non-wetting phase in this work.Hence, we found no noticeable difference in the relative permeability and the relative streaming potential coupling coefficient modelled using either straight or variable radii capillary tubes;  Our model produces considerably different relative permeability curves with small irreducible water saturation (<0.2) in comparison with previously published studies of Jackson [21] and Soldi et al. [37].However, there is no noticeable difference between modelled curves using either of the approaches if irreducible water saturation is larger than 0.2;  Compared with the results published by Jackson [21] and Soldi et al. [37], the relative streaming potential coupling coefficient simulated with our model appears to be more stable at high water saturation and to decrease more rapidly to zero as water saturation approaches the irreducible value.This behaviour is consistent with published experimental results, thus suggesting that the non-monotonic capillary size distribution should be used for more accurate characterisation of multi-phase flow in porous media;  Our model was used to simulate measurements of the streaming potential coupling coefficient in sandstone samples saturated with aqueous solution and liquid CO2 [44].
The model assumed a thick double layer approach for computing the coupling coefficient, consistent with the use of tap water in the experiments.The modelling results suggest that true irreducible water saturation was not reached in the experiments reported by Moore et al. [44].This conclusion is consistent with the hypothesis that explained a non-zero coupling coefficient in the experiments of Vinogradov and Jackson [18].Moreover, since our model produces qualitatively more accurate behaviour of the coupling coefficient with decreasing water saturation, the discrepancy between the reported by Moore et al. [44] irreducible water saturation and the modelled true value was the smallest with our approach relative to that of Jackson [21] or Soldi et al. [37];  To improve the quality of the here-developed bundle of capillary tubes model requires explicit representation of the residual (capillary) trapping of the non-wetting phase.This modification will be developed in a future study using the alternating capillary radii and will potentially allow a more accurate depiction of hysteretic behaviour of the streaming potential coupling coefficient during saturation and desaturation of the modelled rock with the non-wetting phase;  Due to its simplicity, the here-reported and to-be-improved bundle of capillary tubes model can be used to accurately simulate the evolution of the streaming potential coupling coefficient during multi-phase flow in porous media, thus providing an efficient means for a variety of geophysical applications.

Figure 1 .
Figure 1.The bundle of capillary tubes (bold curves) model.The inset to the right shows the two representations of modelled capillaries having the (i) constant (left) and (ii) alternating (right) radius.In the alternating capillary radius model, the pore throat is modelled by a reduced radius , while the pore throat length is defined by , where is the pattern length, which is described in more detail below[35].

Figure 2 .
Figure2.Pore size distribution of Berea sandstone reported in terms of (a) pore radius, extracted from Shi et al.[39], and (b) pore throat radius, extracted from Li and Horne[38].

Figure 5 .
Figure 5. Relative permeability (a-c) and relative streaming potential coupling coefficient (d-f) of the BOCT model constructed using the Soldi (a,d), Jackson (b,e) and New (c,f) distribution functions with capillaries of alternating radius assuming a thick electrical double layer and significant surface conductivity with maximum capillary radius of 100 μm.

Table 1 .
Pore and throat radii of Berea sandstone samples with respective porosity and absolute permeability.NMR is the Nuclear Magnetic Resonance.

Table 2 .
Values required for each parameter of Soldi, Jackson and New CSD functions to obtain porosity of 18.75% for constant radius BOCT model when = 0.1 μ .Number of significant figures varies due to the requirement to match porosity exactly.

Table 3 .
Permeability and number of capillaries within the model for Soldi, Jackson and New CSD functions with constant capillary radius and a model porosity of 18.75% when = 0.1 μm.

Table 4 .
Values required for each parameter of Soldi, Jackson and New distribution functions to obtain porosity of 18.75% with alternating radius water wet bundle of capillary tubes model with radial and length factors of 0.7 and 0.2, respectively, when = 0.1 μm.Number of significant figures varies due to the requirement to match porosity exactly.

Table 5 .
Permeability and number of capillaries within the model for Soldi, Jackson and New distribution functions with alternating capillary radius and model porosity of 18.75% with radial and length factors of 0.7 and 0.2, respectively, when = 0.1 μm.

Table 6 .
[43]ace conductivity of the model for Soldi, Jackson and New distribution functions with capillaries of alternating radius and model porosity of 18.75% required to match the conductivity of the water-saturated rock measured by Moore et al.[43].Number of significant figures varies due to the requirement to match conductivity exactly.

Table 7 .
Actual irreducible water saturation of the Berea sample for various reported irreducible water saturations when ( ) = 0.1, assuming the presence of a thick electrical double layer.Actual irreducible water saturations determined for Soldi, Jackson and New distribution functions with capillaries of alternating radius and model porosity of 18.75%, when = 100 μm.