International Journal of Geophysics
Volume 2019 (2019), Article ID 2985824, 19 pages
https://doi.org/10.1155/2019/2985824

Research Article

Self-Potential Studies in Volcanic Environments: A Cheap and Efficient Method for Multiscale Fluid-Flow Investigations

N. Grobbe1,2 and S. Barde-Cabusson¹

¹University of Hawai‘i at Mānoa, Hawai‘i Institute of Geophysics and Planetology, School of Ocean and Earth Science and Technology, USA
²University of Hawai‘i at Mānoa, Water Resources Research Center, USA

Correspondence should be addressed to N. Grobbe; ngrobbe@hawaii.edu

Received 8 February 2019; Revised 30 June 2019; Accepted 18 July 2019; Published 20 October 2019

Academic Editor: Filippos Vallianatos

Copyright © 2019 N. Grobbe and S. Barde-Cabusson. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

We demonstrate the value of using the self-potential method to study volcanic environments, and particularly fluid flow in those environments. We showcase the fact that self-potential measurements are a highly efficient way to map large areas of volcanic systems under challenging terrain conditions, where other geophysical techniques may be challenging or expensive to deploy. Using case studies of a variety of volcano types, including tuff cones, shield volcanoes, stratovolcanoes, and monogenetic fields, we emphasize the fact that self-potential signals enable us to study fluid flow in volcanic settings on multiple spatial and temporal scales. We categorize the examples into the following three multiscale fluid-flow processes: (1) deep hydrothermal systems, (2) shallow hydrothermal systems, and (3) groundwater. These examples highlight the different hydrological, hydrothermal, and structural inferences that can be made from self-potential signals, such as insight into shallow and deep hydrothermal systems, cooling behavior of lava flows, different hydrogeological domains, upwelling, infiltration, and lateral groundwater and hydrothermal fluid flow paths and velocities, elevation of the groundwater level, crater limits, regional faults, rift zones, incipient collapse limits, structural domains, and buried calderas. The case studies presented in this paper clearly demonstrate that the measured SP signals are a result of the coplay between microscale processes (e.g., electrokinetic, thermoelectric) and macroscale structural and environmental features. We discuss potential challenges and their causes when trying to uniquely interpret self-potential signals. Through integration with different geophysical and geochemical data types such as subsurface electrical resistivity distributions obtained from, e.g., electrical resistivity tomography or magnetotellurics, soil CO₂ flux, and soil temperature, it is demonstrated that the hydrogeological interpretations obtained from SP measurements can be better constrained and/or validated.

1. Introduction

Volcanic environments are fascinating and complex geological settings. Besides activity that is directly visible at or above the surface, a lot of complicated, dynamic processes happen below the surface that determine the subsurface characteristics of the volcano, which in turn control, e.g., groundwater flow and local hydrology [1, 2]. Viewing volcanic settings and their subsurface fluids in light of renewable energy, the following factors are ingredients for a successful geothermal play [3]: sufficient subsurface heat (temperature), the presence of fluid, and high permeability (for fluids to be extracted and replenished).

Volcanic settings are particularly unique in terms of heat anomalies affecting various subsurface fluid-flow processes. For example, both in active and extinct volcanoes, subsurface heat and remnant heat may generate hydrothermal fluid motion [3–5].

Permeability in volcanoes is controlled by several factors. These include the inherent formation properties of the volcanic deposits and possible chemical and/or mechanical alteration (weathering) of formations. The permeability is also heavily influenced by the presence of cracks, faults, and fracture networks. These features can strongly control groundwater flow paths. In volcanic systems, a variety of mechanisms can develop cracks, faults, and fracture networks, which can enhance the permeability and establish preferential fluid flow paths (although faults may also have a sealing-nature). At least four general scenarios seem plausible for developing crack/fault/fracture networks in volcanic environments:

(1) Rift Zones/Dike Intrusions. Although the details of dike formation are a topic of research, a substantial body of research indicates that dikes are controlled by the stress field at the time of intrusion, with dikes propagating along surfaces normal to the direction of least compressive stress [6, 7]. During dike intrusion, the dikes themselves perturb the stress field over distances that scale with the height or length of the dike [7]. Laboratory and field experiments demonstrate that fracturing and other types of inelastic deformation can occur above the top of a dike and presumably elsewhere around its perimeter [8]. The resulting cracks, faults, fractures, and fissures can all enhance the permeability of the volcanic system. The dikes themselves are often considered to be low-permeability barriers for fluid flow [1], although they commonly develop cooling fractures that extend across the dikes; fractures in these intradike fracture systems can be closely spaced [9].(2) Near-Surface Normal Faults and Growth Fault Systems. Basaltic dike intrusions commonly develop in regions of pronounced normal faulting, with networks of nearly vertical opening-mode 1 fissures being commonly along the normal faults [10]. Large voids in these systems can remain open to depths as large as several hundred meters, and they would play an important role in groundwater systems and hydrothermal circulation. Lava flows that encounter these fissures do not necessarily invade all the void space in the fissures, so the lavas can seal the tops of the fissures while allowing substantial interconnected void space to be retained [10].(3) Cooling Cracks. Cooling of basalt is a diffusion process, which causes contraction and crack formation [11]. Relative temperature gradients, and hence cooling rates, are suggested to play an important role in crack growth processes [12]. The features of the cracks seem to be the result of a delicate coplay between thermal and mechanical behavior during cooling. For example, each increment of crack advance can be halted due to, e.g., lowering of the tensile stress associated with thermal contraction due to displacements accompanying the fracture or as a result of the fracture proceeding into hotter, less brittle material [12]. When considering cooling rates of typical intrusions as well as the temperature gradients within intrusions, the environment plays an important role: e.g., groundwater and vapor interactions with the magma [13].(4) Flank Collapse/Landslide-Related Faults. Some volcanic flank collapses are essentially large landslides. The landslide process is well known to develop fractures around the perimeter of a slide mass, as well as in some areas within the slide mass. Sliding at depth can cause fracturing and normal faulting at the surface, as well as thrust faulting near the toe of the landslide [14]. A combination of gravitational forces, magma pressure, and groundwater seepage effects all have been proposed to help drive massive landslides on Hawaiian volcanoes [15]. Phenomena such as hydrothermal convection and hydrothermal alteration also contribute to volcanic slope failure [15].

The main processes in volcanic environments that reduce rock formation permeability (or hydraulic conductivity) are chemical alteration processes [16]. The original hydraulic conductivity of the rocks can be drastically reduced or even totally obliterated through processes such as ferralitic weathering (a soil genetic process that leads to the removal of bases and silica [17]) and thermal and hydrothermal alteration [18, 19]. For example, acidic hydrothermal fluids can chemically alter rock formations, forming, e.g., clay layers with low permeability in certain zones. For instance, strong chemical alteration (weathering) can occur in the vicinity of preferential fluid flow paths and/or nearby magmatic conduits and hydrothermal vents [4]. Chemically altered formations, due to their low hydraulic conductivity, control groundwater flow paths or may form barriers/traps for groundwater (e.g., perched aquifers). Hence, in volcanic environments, a complex coplay exists between permeability enhancing fracture processes and permeability decreasing chemical alteration processes.

The clay minerals, besides reducing the permeability of the rock formation, have a pronounced effect on the electrical conductivity of the medium as well. Due to the relatively high specific surface area of clay minerals, combined with the CEC (cation exchange capacity) of clays, chemically altered formations can possess high surface conductivities [20]. For fully water-saturated, porous clayey materials, we can write the complex electrical conductivity as , where is the electrical conductivity of the pore fluid, is the complex valued surface conductivity, and is the formation factor [20]. Therefore, high electrical conductivities can correspond to chemical alteration fronts [21]. The electrical conductivity of the medium affects the signals of electrical geophysical techniques, such as self-potential measurements, and hence the presence of clays can have a significant impact on measured SP signals.

In the above, we have briefly introduced various levels of complexity that exist in volcanic environments, affecting hydrology: complex geology, a strongly spatially varying and dynamically changing formation permeability/hydraulic conductivity, together with the presence of heat and heat-driven fluid motion. The systems are further complicated by the variability between different types of volcanoes and volcanic features, such as stratovolcanoes, shield volcanoes, tuff cones, and monogenetic volcanoes.

These features make volcanic environments very interesting and challenging study sites, especially in terms of groundwater hydrology and geothermal energy.

In terms of (hydro)geophysics, volcanic environments are highly challenging environments to image and characterize. For example, basaltic rocks typically result in very high seismic velocities, and seismic imaging suffers from high amounts of scattering due to the abundance of fractures, cracks, and faults [22].

Ambient noise surface wave tomography has been proven a valuable tool to identify structural subsurface characteristics of volcanic systems, both laterally and with depth [23, 24]. Furthermore, in recent years, the use of ambient noise for monitoring of subsurface fluids is actively being explored [25]. The cost of operation for ambient noise studies is much lower than that for active seismic surveys, since no active, human-made sources are to be deployed. For ambient noise tomography, the distribution of noise sources (affecting illumination) is an important factor, whereby uneven distributions can introduce a bias in the phase and/or group velocities [26].

Electrical Resistivity Tomography (DC resistivity) can provide good quality results for studying volcanic environments [4]. The method actively injects current into the ground using two current electrodes and measures the resulting potential differences at other pairs of electrodes. In this way, an image of the distribution of electrical resistivity in the subsurface can be obtained. However, hydrogeologic interpretation difficulties may arise when trying to distinguish fluid-saturated zones from chemical alteration (e.g., clay layers), both displaying relatively high electrical conductivity. One way to try and distinguish the conductivity associated with the pore fluid from the conductivity related to chemical alteration is via induced polarization, where variations in chargeability of the rock formations are being mapped [4].

The image resolution that can be obtained through static ERT imaging is often insufficient to distinguish individual fractures [27]. Additionally, its depth of investigation can be limited depending on various factors, such as subsurface electrical conductivity distribution, type of equipment and acquisition geometry, and surface contact resistance (difficulty to inject the current) (e.g., [28]). In terms of operation, ERT surveys can be quite time consuming, especially when large depths of investigation are desired (i.e., long cables, large current electrode separations). With cabled, 2D systems, 2D and pseudo-3D resistivity images can be obtained. In recent years, novel equipment has been developed, such as the FullWaver system (http://www.iris-instruments.com/), which is a nodal-based system that measures both the horizontal components of the electric field, thereby allowing for true 3D images. The system can operate in different acquisition configurations, including pole-dipole and dipole-dipole. The pole-dipole capability of the system combined with externally powered current transmitters allows for greater depths of investigation.

An alternative to ERT is the Time-Domain Electromagnetic (TDEM) method (or transient electromagnetics TEM). Electrical currents are induced in the Earth via electromagnetic induction. TDEM often offers better depth resolution than ERT and often displays a higher sensitivity to conductive layers [28, 29]. However, it is more difficult for TDEM to resolve the resistivity of very resistive zones above a conductor [28]. The method allows for faster acquisition over larger areas, avoids challenges with current injection into a resistive surface layer, but is typically more expensive to operate.

In volcanic systems, groundwater and geothermal resources can reside at depths > 2 km (e.g., 2800 meters below sea level [30]). Particularly for volcanic island settings, which are dependent on their own fresh water resources, these deeper sources of groundwater and geothermal energy may be valuable future resources.

The magnetotelluric (MT) method is one of the methods, besides seismics, that can be used to study these deep fluid systems. MT allows for great depths of investigation, but typically at the cost of resolution [31–33]. MT is a passive method which measures both electrical and magnetic fields related to naturally occurring currents. It can cover large areas, in both 2D and 3D. However, the time and effort to deploy a single station is quite significant, especially in volcanic environments where the necessary digging and burial of the measurement instruments is extremely challenging.

In recent years, we have seen the emergence of airborne geophysical methods, such as airborne electromagnetics (e.g., [34–37]). These methods can offer tremendous areal coverage, but typically have very high operational costs as well.

Throughout the years, the self-potential (SP) technique has been proven to be a cheap, efficient, and quite successful geophysical method for studying subsurface fluid behavior in volcanic environments (e.g., [38–42]). As a consequence of the electrokinetic mechanism that generates SP signals (see more below), the SP method can provide unique insight into the actual fluid flow-paths [43] and distribution of Darcy flow velocities [44], thereby assisting in the identification of fracture networks and other zones of high or low permeability/hydraulic conductivity. The method allows for great acquisition flexibility at high speed, allowing for coverage of large areas of land, and the equipment is very cheap.

In this paper, we present an overview of a variety of self-potential field studies carried out in different volcanic environments, including shield volcanoes, monogenetic volcanoes, tuff cones, and stratovolcanoes. Based on these different case studies and the interpretation of the measured SP signals, we illustrate and discuss the variety of volcanic subsurface knowledge, as well as the breadth of information on porous medium fluid processes that we can extract and infer from these datasets, as well as highlight the demonstrated successes of using the SP method for studying volcanic environments. Like other geophysical methods applied to volcanic systems, the SP method has its own benefits and drawbacks. We will discuss several of these in greater detail below, thereby also discussing some of the challenges and pitfalls for SP in volcanic environments.

2. Method

2.1. Brief Background on SP Field Data Acquisition and Processing

The self-potential method measures the naturally occurring difference in electrical potential between two electrodes coupled to the Earth. It is a passive geophysical method; i.e., it does not use human-made sources for the geophysical data acquisition. Different types of nonpolarizable electrodes have been and are currently used, consisting of a metal rod immersed in a solution of this metal. Currently, the most frequently used types are the Pb-PbCl₂ [45, 46] and the Cu/CuSO₄ electrodes. The difference in electrical potential is measured by connecting the electrodes to a voltmeter through an insulated wire of a few hundred meters length. The voltmeter itself is a calibrated, high-impedance voltmeter (approximately 10 to 100 MOhm) with a sensitivity of 0.1 mV. The so-called reference electrode is placed at the beginning of the profile, and each measurement point is performed by digging a small hole, where the moving electrode measures the potential difference between itself and the reference electrode. Prior to each SP measurement, it is recommended to check the resistance between the two electrodes through the ground. This allows for quality control of the electric contact, before measuring the electric potential difference corresponding to the SP value. An alternative technique is to measure the potential difference several times over a small area, rejecting strongly anomalous values within those measurements and assigning the average value to the measurement point. Using the leapfrog method, the reference and moving electrodes exchange roles every few hundred meters in order to cancel out a systematic error build-up due to a possible difference in electrical potential between the electrodes themselves (not through the ground), which is an artifact. Data acquisition is carried out by forming closed loops that create a completely interconnected network [20]. With respect to data processing, the whole dataset is set to a unique reference, in particular by applying a closure correction along these loops, to limit cumulative errors caused by, e.g., changes of environmental parameters or a drift of the SP measurement devices.

The self-potential method offers very fast and low-cost data acquisition, enabling surveying of large areas of land [47]. It is a highly mobile, light-weight technique, allowing for data acquisition in rough and challenging terrain, such as volcanoes.

2.2. Brief Theoretical Background for SP: Microscale Origin vs. Macroscale Processes

Several transport processes can occur in porous media, including coupled processes, i.e., flow as the result of an actuating gradient of a different type than the flow phenomenon itself [48].

Onsager [49, 50] discovered that there is an underlying symmetry for the phenomenological coefficients between actuating force and resulting flux, also known as the Onsager reciprocal relations. Table 1 displays the so-called Onsager matrix, which shows the types of coupling that can occur in a porous medium containing a fluid.

Table 1: Onsager matrix displaying the types of coupled flow that can occur in a porous medium. The main diagonal represents direct flow phenomena, i.e., the actuating force and the corresponding flow are of the same type. The off-diagonal entries represent the coupled flow phenomena. Onsager reciprocity is reflected by the symmetrical entries across the main diagonal (e.g., streaming potential vs. electroosmosis). Bear in mind that this table is a highly simplified representation. Table 1 is modified from [48].

Flux
Driving force: gradient of the actuating phenomenon

Hydraulic
Electric
Chemical
Thermal

Fluid
Darcy’s law
Electroosmosis
Chemical osmosis
Thermoosmosis

Charge
Streaming potential
Ohm’s law
Dorn potential
Seebeck effect

Solute
Ultrafiltration
Electrophoresis
Fick’s law
Soret effect

Heat
Isothermal heat
Peltier effect
Dufour effect
Fourier’s law

Self-potential signals can arise from a variety of these coupled flow phenomena.

Often, the electrical double layer (EDL) that exists at the surface between the solid particles (e.g., grain surfaces in a rock formation) and the pore water [20] plays an important role in the generation of self-potential signals. In a porous medium containing a fluid, a net charge is created on the mineral surface as a consequence of chemical reactions between the mineral surface and the pore water. This surface charge is compensated by the combined effect of the Stern layer [51] and the diffuse layer [52, 53], which together form the EDL. The Stern layer consists of adsorbed ions on the mineral surface [51], whereas the ions of the diffuse layer obey Boltzmann statistics in the Coulombic field that is created by the deficiency or excess of electrical charges on the mineral surface and the Stern layer [52–54].

When the grains and the pore water are moving, the place of zero relative velocity is called the shear-plane. This shear plane occurs somewhere in the diffuse layer. For the sake of simplicity, it is often considered to be at the interface between the Stern and diffuse layer. The electrostatic potential on this plane is referred to as the zeta potential (see, e.g., [20]).

As a consequence of the presence of an EDL, electrokinetic phenomena can occur, such as the electroosmotic phenomenon, where an electric field drives fluid flow, and the streaming potential phenomenon, where a hydraulic gradient creates an electric field (see Table 1).

Various electrochemical effects can induce SP signals. A common mechanism involves the coupling between chemical gradients and electric current density (diffusion of ions, Fick’s law) [55, 56]. Another mechanism is associated with redox processes, e.g., related to ore bodies (mineral exploration) and contaminant plumes [55, 56].

Besides streaming potentials and electrochemical effects, thermoelectric effects (such as the Seebeck effect) can occur, where direct conversion from temperature gradients to electric voltages takes place. The electrical double layer has its own contribution to the total thermoelectric coupling [57, 58]. Additionally, thermoelectric effects can distort the SP measurements [57, 59].

When we interpret SP signals, we often assume that they are predominantly the result of the electrokinetic streaming potential phenomenon [60]. However, given the dynamic nature and complexity of volcanic environments, it is important to be aware of the various physical mechanisms at the microscale that can cause, and influence, self-potential signals and measurements.

When only considering SP signals of electrokinetic origin, different macroscale mechanisms can be responsible for the measured signals, for example, topographic effects and hydrothermal circulation effects. Meteoric water originates in the atmosphere and reaches the Earth’s saturated zone through infiltration and percolation processes [61]. Due to the wide variety in hydraulic conductivity in volcanoes, there exist variations in rainfall infiltration to different geological beds, whose topography governs gravitational flow [61]. In a simplified scheme, the self-potential value generally increases with decreasing elevation, referred to as the topographic effect. This effect is a consequence of the distance along which the downward fluid flow travels from the topographic surface to the saturated zone; that distance typically increases proportionally with elevation [38, 62].

In active volcanoes, shallow thermal energy can trigger hydrothermal circulation that will influence the topographic effect [61]. The combination of upward flow above warm thermal anomalies (e.g., magma reservoirs) and along tectonic interfaces (such as crater walls, faults) due to thermal sources, resulting in a positive SP anomaly, and downward flows associated with counter flows and infiltration processes, resulting in negative anomalies [61], is considered a plausible explanation for the typical W-shape SP profiles observed on active volcanoes [41].

When measured in the field, especially in volcanic environments, SP signals are likely the result of a mixture of various microscopic and macroscopic mechanisms, where the relative importance of individual mechanisms varies for different scenarios. Think for example of SP anomalies related to geothermal systems: strong thermal gradients exist in the subsurface, and much of the geothermal energy is released through groundwater circulation, hot gas emission, and thermal conduction [61]. Even though it may be challenging to uniquely identify the individual contributions of each mechanism to the measured SP signals, the fact that SP signals are sensitive to both the thermal state and hydrothermal state of volcanoes makes them a unique tool for studying these volcano dynamics [61].

The above is to emphasize that besides the variety and/or combination of microscale physical processes that can generate SP signals (e.g., electrokinetic, electrochemical, and thermoelectric), the actual macroscale structural features play a major role in the characteristics of the measured SP signals. Measured SP signals are the result of the coplay between microscale processes and macroscale structural and environmental features.

For the volcanic case studies presented in this paper, the measured SP signals are primarily interpreted as being a consequence of macroscale features and flow behavior and assume the electrokinetic (streaming potential) mechanism to be the major microscopic origin for these signals.

3. Case studies

We present several self-potential case studies and their accompanying hydrogeological interpretations for a variety of volcanic environments, including tuff cones, shield volcanoes, stratovolcanoes, and monogenetic fields.

These studies showcase the success and value of SP data acquisition for studying volcano hydrogeology and dynamics and highlight the applicability of the method for these various volcanic environments.

We categorize the examples into the following three multiscale fluid-flow processes:

(1)Deep hydrothermal systems: La Fossa, Piton de la Fournaise, and Teide(2)Shallow hydrothermal systems: Piton de la Fournaise and Stromboli(3)Groundwater: Teide and Garrotxa.

These examples display SP signals associated with several of the mechanisms discussed above, including

(i)Volcanic hydrothermal systems displaying the “W”-shape SP signal pattern(ii)Shallow hydrothermal anomalies, from the discharge of the central hydrothermal system and from recent volcanic deposits(iii)Thermoelectric effects related to the cooling of lava flows(iv)Structural and hydrogeological features: crater limits, regional faults, rift zones, incipient collapse limits, different hydrogeological domains, elevation of the water level/aquifer, and buried caldera.

3.1. Deep Hydrothermal Systems

3.1.1. Tuff Cone (Stratovolcano): La Fossa Di Vulcano, Vulcano, Aeolian Islands, Italy

Our first example is the stratovolcano La Fossa (Vulcano, Aeolian Islands, Italy). The volcanic activity of Vulcano is caused by magmatic activity that is controlled by regional tectonics, such as strike-slip faults [63]. La Fossa is a relatively small edifice and has a strong hydrothermal character [64], and it is located where the most recent volcanic activity of the island took place. North of the Fossa cone is the small Vulcanello edifice, which forms a roughly circular peninsula that connected to the main island in 1550 AD, during its last eruption. The latest eruption from Vulcano consisted of explosive activity from the Fossa cone from 1888 to 1890 (e.g., [65]).

Active for ~6000 years [66, 67], La Fossa cone is the current active volcanic center of the Island of Vulcano. Built inside the Fossa caldera, which was formed after the collapse of the Cardo tuff cone and the Lentia intrusive Complex, the present cone is the result of several phases of activity [68]. Surges, fallout deposits, and effusive activity built up the Fossa cone with stacked and intersected ‘a‘ā and rhyolitic flows, intrusions, pyroclastic ashes, and tuff, which are all affected to different degrees by an intense hydrothermal activity and by regional tectonics, as well as by at least five major craters associated with various explosive eruptions. This makes La Fossa di Vulcano a complex and highly heterogenous volcanic edifice.

Figure 1 (top) displays the SP data (black dots) as collected along a 2D transect over the cone [64]. Figure 1 (bottom) shows the result of a 2D electrical resistivity tomography (ERT) profile [64]. We can clearly observe a positive SP anomaly coinciding with the electrically conductive zone in the center of the crater (blue colors). The low electrical resistivity values in this central part of the cone are most likely caused by the presence of a hot, hydrothermal system combined with chemical alteration of the deposits (resulting in various clay minerals).

Figure 1: 2D measured and modeled self-potential profile (top) and 2D electrical resistivity tomography transect with modeled groundwater flow velocities (bottom), La Fossa di Vulcano, Vulcano Island, Italy. We can clearly observe the positive SP anomaly coinciding with the electrically conductive zone in the center of the crater. Finite-element modeling of the groundwater flow using Comsol Multiphysics 3.3 confirmed that the observed SP anomalies are very likely to be caused by the combination of a hot, hydrothermal system in the center, causing hydrothermal upwelling and a positive SP anomaly, and downward flow through ashes and tuff materials along the flanks of the volcano. The figure is modified from [64].

As a consequence of hydrothermal upwelling, a positive SP anomaly can be observed (Figure 1, top, red-shaded area).

On the flanks of the cone, the negative SP anomalies are interpreted to be the result of downward flow and infiltration/percolation of cooling groundwater at relatively shallow levels (groundwater flow paths are controlled by structural and geological features of the volcano). Both interpretations have been validated by finite element modeling of the groundwater flow paths, using Comsol Multiphysics 3.3 [64]. The arrows in Figure 1 (bottom) represent normalized groundwater/seepage velocities and highlight the directionality of the groundwater flow, which matches the upwelling and downward percolation interpretations [64]. Based on these groundwater flow paths, the corresponding SP anomalies were calculated, displayed by the black solid line in Figure 1 (top). It highlights that the groundwater flow patterns as modeled indeed explain the observed SP anomalies over the La Fossa cone [64]. One can observe that the two red, resistive anomalies act more or less as no-flow bodies. This is a reasonable result, highlighted by the good match between modeled and observed SP values. Moreover, the high electrical resistivity values of the red anomalies match those of dry basalt (e.g., [31]), thereby indicating that very little fluid resides/flows through these structures.

Figure 2 shows the SP map of the whole cone and illustrates the complexity faced when interpreting only SP data. Here, some strong positive and negative anomalies visible in the SP map remain uncorrelated with other data types collected in the area (i.e., soil temperature and soil CO₂ flux). In the case of Vulcano, the corresponding datasets were acquired at three different times during a one-year period [63]. Volcanic activity and atmospheric conditions vary significantly year-long. Also, sharp lithological contrasts inside the edifice and heterogeneities related to the hydrothermal alteration of volcanic deposits may affect the zeta potential and consequently the measured SP signal in a heterogeneous manner, producing artifacts in the final SP map. In such a scenario, the most reliable information is then obtained in areas with steady hydrothermal activity. This is for example the case in the central crater where the main hydrothermal system has its major surface expression, as well as in former craters, along crater faults, and at locations of diffuse activity inside those craters. In such a scenario, it seems more convenient to rely on SP profiles (Figure 3), which guide interpretations at a smaller scale and can be more easily correlated to data types obtained from other methods at the same location, such as a 2D ERT profile. In La Fossa di Vulcano, this approach allows to define structural limits, such as the limits of buried craters affecting subsurface fluid flow. In any case, interpretation of SP data always benefits from a joint study with results obtained from other data types. This can be for example basic knowledge of the geological setting of the study area, joint temperature or soil gas-detection measurements, or 2D/3D geophysical imagery (e.g., ERT, electromagnetic methods, radar, etc.) and can provide additional insight into the source and the pathways of the fluids responsible for the SP signal. The different methods and data types will inform each other, thereby refining final interpretations.

Figure 2: Self-potential (top) and soil temperature (bottom) maps of La Fossa cone overlain on the Digital Elevation Map (DEM) of the area. White dots represent the measurement points. Areas with steady hydrothermal activity are well delineated by self-potential maxima and correlated to high temperature anomalies inside the craters and along their limits. PN (Punte Nere), Pa (Palizzi), FV (Forgia Vecchia), PC (Pietre Cotte), and GC (Gran Cratere) craters are localized with black-dashed lines. Profile marked with two orange dots corresponds to the cross sections and transects of Figures 1 and 3. This figure has been modified from [63].

Figure 3: Soil temperature measurements (top), self-potential values (middle), and electric resistivity tomography results (bottom) along Profile (see Figure 2). Compared to SP maps (of for example Figure 2), the data visualization in terms of single profile-values allows for a finer-scale analysis of the signal and a more meaningful comparison with other data types. PN: Punte Nere crater, PC: Pietre Cotte crater, and GC: Gran Cratere crater. This figure has been modified from [63].

The Fossa cone example demonstrates that SP measurements are capable to (1) detect the presence of hydrothermal systems on active volcanoes, (2) identify the percolation and infiltration of groundwater along the flanks and in shallow regions of volcanic edifices. Similar observations have been made in Stromboli [5], and (3) inform on subsurface fluid flow paths, especially when combined with ERT surveys and groundwater modeling.

3.1.2. Shield Volcano/Volcanic Island: Piton De La Fournaise, France

The Piton de la Fournaise volcano is located on La Reunion Island, in the Indian Ocean. With 2 or more eruptions per year since 1998, it is one of the most active volcanoes in the world. Most of the eruptive activity takes place in the depression formed by the Enclos Fouqué caldera and the Grand Brûlé. This U-shaped depression is open to the East (Indian Ocean) (Figure 6) and it is thought to be the result of a collapse, flank slide, or a combination of both processes. The Bory crater, one of the two summit craters, did not experience major volcanic or tectonic activity of its own [69]. Looking only into the recent history of Piton de la Fournaise, the Bory crater remains almost unchanged since at least 1766, and only experienced a few minor eruptions (in 1937, 1942, 1956, 1979, and 1981), fracturing, and the collapse of its most eastern wall in 1953 [69]. On the other hand, the Dolomieu crater has already experienced numerous collapse events and the crater floor was repeatedly covered by new lava flows and projections [69]. In April 2007, the Dolomieu crater was finally affected by a major collapse [70, 71] and subsequent eruptions are now slowly refilling the crater.

On Piton de la Fournaise, the whole terminal cone was mapped using the SP method, with a measurement spacing of 20 m [72] along the profiles, creating a dense dataset, and the measurement profiles were further densified in the Dolomieu crater (Figure 4). The SP mapping of the terminal cone highlights different areas displaying high SP values that are located along structural limits of the volcano. Those structural features have thus been interpreted as preferential flow paths for uprising hydrothermal fluids. The most striking examples here are the high SP values highlighting the well-known rift zones of Piton de la Fournaise, which are intersected by a dense fault-network. The summit Bory and Dolomieu craters clearly show SP maxima (yellow to red colors) as well and act as barriers for hot fluid circulation. When focusing on the summit craters (Figure 4), smaller structures such as Pre-Bory, Soufrière, and Petit Plateau, former pit craters surrounding the Dolomieu, are also highlighted by SP maxima. The authors show as well how SP mapping reveals fluid flow spatial heterogeneities inside the Dolomieu crater. As the Dolomieu collapsed during the April 2007 eruptive crisis of Piton de la Fournaise, this heterogeneity eventually proved to be the expression of a fault controlling the collapse sequence: the crater floor section with lower SP values (due to less hydrothermal activity) was temporarily preserved during an early stage of the collapse, forming an eastern plateau [72].

Figure 4: Self-potential map of the whole terminal cone of Piton de la Fournaise (top) and a detailed SP map of the crater (bottom). PB (Pre-Bory), S (Soufrière), and PP (Petit Plateau) are former pit craters surrounding the Dolomieu and Bory craters, affected by major hydrothermal release as shown by the high self-potential values. White-dotted lines represent the collected SP profiles. Self-potential maxima highlight the rift zones of this volcano in consistency with the repartitioning of the fractures and fissures. Before the 2007 eruption, the self-potential data identified underground fluid flow heterogeneities highlighting structures that later controlled the collapse of the Dolomieu crater and the temporary preservation of an eastern plateau. Figure modified from [72].

In the case of Piton de la Fournaise, SP has shown to be a powerful tool to highlight hydrothermal circulation heterogeneities controlled by rift zone fracturing processes, crater limits, and incipient collapse faulting.

3.1.3. Stratovolcano/Volcanic Island: Teide, Tenerife, Canary Islands, Spain

Our next example is the Central Tenerife Volcanic Complex (Canary Islands, Spain). The central area of Tenerife is dominated by Las Cañadas caldera where two stratovolcanoes, Pico Viejo and Teide, are located. The initial stages of activity of this island resulted in the formation of a basaltic shield which emerged from the ocean more than 11.9 Ma ago [73]. The Las Cañadas edifice was built with more evolved magmas (mostly basalts, trachytes, and phonolites) and experienced several cycles of construction and destruction between 4 Ma and ~0.17 Ma, while more basaltic activity occurred along a NW rift zone, a SE rift zone, as well as a diffuse southern basaltic field [74].

The Las Cañadas edifice eventually collapsed forming Las Cañadas caldera at 0.18-0.13 Ma [75]. It is composed of three smaller calderas (Ucanca, Guajara, and Diego Hernández) and various lateral landslides were identified in the area. A debate exists concerning the mechanisms involved in the formation of this 16 km wide elliptical caldera, whether due to vertical or horizontal offset or a combination of both [76]. Pico Viejo (3135 m a.s.l.) and Teide (3718 m a.s.l.) further built up in the caldera with emission of basanites, tephrites, and phonolites, as well as occasionally more evolved magmas [77]. During the last 180 ka, mafic to intermediate eruptions also occurred in the central zone, along the rift zones, and with minor activity in the southern basaltic field as well [74, 78].

Beyond the fumarolic activity in some parts of the caldera floor and on the two stratovolcanoes, the most recent activity of the Tenerife volcanic complex was a seismic crisis in 2004, which consisted of a swarm of more than 1000 low-magnitude earthquakes of magmatic origin [79]. However, no eruption accompanied this crisis.

For studying the volcano and its fluid dynamics over such a complex and huge area, the SP method is probably the most suitable land-based geophysical tool, as it offers the possibility to cover long distances with high resolution and in a small amount of time compared to other methods. Figure 5 presents the results of an extensive SP survey carried out in Las Cañadas caldera. It is shown that the SP signals allow to identify and map the Teide and Pico Viejo hydrothermal systems and to delineate areas with different hydrothermal intensity. The zones of different hydrothermal intensity at these two stratovolcanoes are thought to be correlated with unknown structures (e.g., former craters, calderas) inside the edifices. The SP maxima also highlight smaller hydrothermal systems at two former volcanic vents (Montaña Blanca and Montaña Rajada). Those may be either active or fossil hydrothermal systems according to [41], who showed that fossil hydrothermally altered volcano cores can display weak positive SP anomalies under certain conditions.

Figure 5: SP survey in Las Cañadas caldera (Tenerife, Spain). The SP maxima highlight the main hydrothermal systems of Teide and Pico Viejo volcanoes, as well as smaller hydrothermal systems at two former volcanic vents (Montaña Blanca and Montaña Rajada). Ground water infiltration is highlighted by self-potential minima at the Diego Hernandez (DH) buried caldera and at the NW rift fractures (NWR). Figure from [76].

Figure 6: Plot of the time-lapse self-potential measurements performed at the coastline of Piton de la Fournaise. The upper graph shows self-potential data acquired in 2004 (green curve) and in 2006 (black curve) as well as the difference between both surveys (bottom black line). Recent lava flows are represented by colored bands located in the DEM figure (bottom right) adopting the same color coding. The figure is modified from [80].

3.2. Shallow Hydrothermal Systems

Despite the potential impact on deep-source interpretations of SP values, shallow sources of SP signals are definitely of interest when studying hydrothermal circulations in active volcanoes. Both deep and shallow hydrothermal systems can have an impact on the measured SP signals. This can occur simultaneously, or, separately, depending on the physical target of interest and the location of the survey on the volcano.

3.2.1. Shield Volcano/Volcanic Island: Piton De La Fournaise, France

Let us revisit the Piton de la Fournaise example that we have discussed previously in terms of deep hydrothermal systems. We now look at time-lapse SP data acquired on Piton de la Fournaise.

Figure 6 shows SP profiles acquired at two different times (i.e., one in February 2004 and the second in July/August 2006) along the same track, in the coastal area of Piton de la Fournaise [80]. Furthermore, it shows the difference in SP values between both surveys. These profiles were acquired along the RN2 national road in the Grand Brûlé, an area frequently affected by lava flows. The plot shows that the two SP signals are similar except for areas affected by recent lava flows. Measurements performed on lava flows crossing the road before February 2004 show a decrease of the SP signal, while those performed on lava flows crossing it between February 2004 and July 2006 show an increase of SP. In this case, SP measurements are able to detect very shallow hydrothermal activity inside young lava flows (of a few decades old). The streaming current associated with the flow of pore water is likely responsible for the high SP values measured on those recent lava flows, and the thermoelectric effect generating SP signals in addition to the electrokinetic streaming potential probably contributes significantly to the signal strength as well [80]. In turn, the decrease in temperature during the cooling of the lava flows, and consequently a decrease in shallow hydrothermal activity inside the flow after their emplacement as well as a decrease in thermoelectric contributions, will result in a decrease of the SP values. This opens the door for a variety of SP measurement applications, ranging from mapping hydrothermal activity to monitoring of hydrothermal activity and lava flow cooling processes over time. On the flip-side, it also means that, in general, determining the source of the measured SP signals uniquely can be challenging (especially when SP data is the only data collected), as the contributions from a shallow SP source will stack onto, and interfere with, contributions from a deeper source. The result will be that the shallow source contribution acts as a parasitic signal in studies aiming to map deep hydrothermal activity and structures, with the risk of it causing misinterpretations [80]. Combining SP data with other geophysical and/or geochemical data types, or at the very minimum with basic knowledge of the geology of the study area, might thus be indispensable in cases of rapidly varying hydrothermal activity, shallow heat sources, and other shallow sources of SP signals.

3.2.2. Stratovolcano/Volcanic Island: Stromboli, Italy

The second example of shallow hydrothermal activity highlighted by SP data is the Stromboli Volcano, in Italy.

Stromboli stratovolcano is the northernmost island of the Aeolian archipelago, located in the Tyrrhenian Sea. The subaerial volcanic history of Stromboli is subdivided into seven main phases, including Paleostromboli I, II, III, Scari complex, Vancori, Neostromboli, and Recent Stromboli (e.g., [81]). Those phases are separated by quiet erosional stages or by caldera or flank collapses (e.g., [82]). The hydrothermal activity of this volcano has been well described by [40], based on SP and CO₂ soil concentration measurements. The authors describe a complex hydrothermal system mainly centered at the summit area but with a minor subsystem slightly offset to the main one.

Looking at the current eruption dynamics at the summit vents of Stromboli, its activity has been almost continuous since 1932. The ordinary activity consists of frequent Strombolian explosions (approximately every fifteen minutes to every hour), as well as occasional lava flows. Most of this activity occurs in the summit area and in a restricted collapsed zone named Sciara del Fuoco, but occasional stronger explosions or paroxysms [83] can throw volcanic bombs up to a few hundred meters to a few kilometers away from the summit vents. One of those events took place on March 1 2007, ejecting one of the volcanic projectiles from the summit active vents (i.e., at ~750 m a.s.l.) all the way to the flank of the volcano at 560 m a.s.l. [5]. After this event, the area of impact began to emit vapors and rose in temperature even a long time after the impact, which is unusual for a volcanic bomb impact crater. The heat from a volcanic bomb usually lasts a few hours or days, depending on the size of the projectile and then dissipates. This particular impact triggered a new persistent vapor emission at this distal location of the Stromboli volcano. This is why this particular event was studied. Figure 7 shows the results of an investigation around the block impact crater in May 2008, using SP, soil temperature, soil CO₂ flux, and electric resistivity tomography. At the location of the impact, the plots display a clear maximum value of SP, correlated to both the temperature and CO₂ flux signals. The resistivity model shows a shallow conductive (<200 Ω.m) layer at 10-15 m depth, nearly parallel to the topographic surface, and connected to the surface through a narrow conductive channel at the location of the anomalies detected by SP and the other methods. This conductive layer was interpreted as shallow hydrothermal activity fed both by (1) the central hydrothermal system of Stromboli, as water condensing in the summit area of the volcano flows downward through more porous, shallow stratigraphical layers, and (2) fluids rising from a deeper source along a regional fault [5].

Figure 7: Soil temperature, self-potential, soil CO2 flux measurements, and electric resistivity tomography results of the multidata study carried out on the Stromboli volcano (Italy). At the location of a volcanic bomb impact crater from the March 15th 2007 paroxysmal eruption, we observe a clear maximum value for the self-potential, which correlates with the other measurements. This study uncovered shallow hydrothermal circulation which seems a common feature of active volcanoes. Figure from [5].

The contribution of the SP method is crucial here as the conductive layer revealed by the electric resistivity tomography alone could alternatively be interpreted as a chemically altered layer rich in conductive clay minerals.

3.3. Groundwater

Thus far, we have looked at the use of the self-potential method to study deep and shallow hydrothermal systems and its fluid circulation in volcanic environments. We will now provide two examples where SP signals inform us on groundwater flow in volcanic systems and the associated hydrogeology.

3.3.1. Stratovolcano/Volcanic Island: Teide, Tenerife, and Canary Islands, Spain

The Teide volcano, previously discussed to refer to its deep hydrothermal system based on an extensive SP dataset published by [76], also shows interesting results from a hydrogeological point of view. The SP map displays hydrogeological signals interpreted as groundwater infiltration at the Diego Hernandez buried caldera and at the NW rift fractures (small SP signals in bluish colors). This is another useful application of the collected SP data. Figure 8 shows the results of the conversion of SP data to a model of the topography of the water table in Diego Hernandez [76]. Previously published audiomagnetotellurics data [84] together with water-level data measured at two wells in the same area were used to convert the SP map into a map of water table elevation. It seems that the underground topography of the caldera floor is a major controlling factor of the hydrogeology in this area. The SP signal demonstrates as well the connection of the hydrogeological systems of the Diego Hernandez caldera and the Orotava valley.

Figure 8: Conversion of the self-potential data (using previously acquired audio-magnetotelluric and water table data) into information on the elevation of the water table. Different colors represent different elevations of the bottom of the vadose zone in meters above sea level. Self-potential demonstrates that the underground topography of the Diego Hernandez caldera floor is a major controlling factor for the hydrogeology in this area. Figure from [76].

3.3.2. Monogenetic Volcano: Garrotxa, Spain

Our final study site is located in Catalonia (Northern Spain). This case is different from the rest of the examples presented in this paper for two main reasons: (1) this is a monogenetic volcanic field that consists of over 50 monogenetic volcanoes and (2) the last volcanic eruption (Rocanegra eruption) is known to be less than 11.5–13 ka old [85, 86]. Monogenetic volcanism is a peculiar type of volcanism as its first characteristic is to have no central path for the magma to reach the surface. Each new eruption arises at a different location, builds an independent volcanic edifice, and lasts from a few hours to a few years. In the Garrotxa, despite a relative monotony in the geochemistry of the erupted lavas, an important variety of eruptive dynamics took place, from lava flows to violent phreatomagmatic explosions. Even during a single eruption, the eruptive style can change, leading to the formation of complex edifices. This monogenetic field lies on top of Upper Palaeozoic granites and schists or on Eocene and Quaternary sedimentary subtracts. Eruptive products are a suite of intracontinental leucite, basanites, nepheline basanites, and alkali olivine basalts [87]. Regional tectonics play an important role on the spatial repartitioning of the volcanism in the field. Eruptive fissures are usually associated with NNW–SSE and NNE–SSW Neogene tectonic structures [88].

In this case study, SP is expected to provide information on the hydrogeological system of the area, structural information, and their relation with volcanism. Figure 9 presents both the obtained SP map of part of the Garrotxa volcanic field [88] as well as an analysis of the measured SP values versus elevation [89]. The SP map shows a significant difference in self-potential values between the western and eastern zones, which suggest a lateral heterogeneity of the medium. The western and eastern zones are separated by a NNW–SSE lineament and a NNE–SSW lineament which could at first be interpreted as the eruptive fissure associated with Santa Margarida and Croscat monogenetic volcanoes. However, Figure 9 (right) drives to another interpretation. In purely hydrogeological zones, SP values are strongly controlled by the vertical distance between the measurement point and the water table, i.e., the local thickness of the vadose zone. Typically, there is a negative correlation between SP and elevation; considering a flat continuous aquifer and a homogeneous medium, SP values decrease with elevation as the thickness of the vadose zone increases (the topographic effect). In the Garrotxa, the eastern zone shows a linear trend with an SP versus elevation gradient of -0.19 mV/m, typical for hydrogeological zones on volcanoes [90]. In the western zone, no particular trend is observed. This would favor the hypothesis of the two areas having distinct hydrogeological properties and/or different global electrical resistivity values, separated by the NNW–SSE and NNE–SSW lineaments.

Figure 9: Self-potential map of the Garrotxa monogenetic field (left). The map highlights the strong lateral heterogeneity of the medium resulting in a complex setting. A western zone and an eastern zone appear to be separated by a NNW–SSE lineament and a NNE–SSW lineament. A plot of self-potential vs. elevation (right) seems to confirm this separation by displaying a strong difference in SP behaviour between the east and west sides, which is likely due to lithological and hydrogeological contrasts between the two areas. Figure modified from [89].

4. Discussion and Conclusions

We have shown that SP measurements are a highly efficient way to map large areas of volcanic systems under challenging terrain conditions, where other geophysical techniques may be challenging or expensive to deploy. Given the complex dynamics of fluid flow, heat, and permeability variations in volcanic environments, the sensitivity of SP data to fluid flow provides a unique way to identify and highlight infiltration patterns, heat-induced shallow and deep hydrothermal convection, and preferential flow paths such as fractures, faults, and impermeable zones.

The various case studies presented in this paper demonstrate that SP signals, despite the fact that they can result from a variety of microscale mechanisms and macroscopic features and consequently can be challenging to uniquely interpret as a stand-alone data type, provide valuable insight into the local hydrogeology, fluid flow dynamics, and macroscale structural features of volcanic systems, regardless of their type, and on a variety of spatial and temporal scales, such as information on

(1)Deep hydrothermal systems(2)Shallow hydrothermal anomalies, from the discharge of the central hydrothermal system and from recent volcanic deposits(3)Cooling behavior of lava flows(4)Different hydrogeological domains(5)Upwelling, infiltration, and lateral groundwater and hydrothermal fluid flow(6)Elevation of the water level/identification of aquifers(7)Crater limits, regional faults, rift zones, incipient collapse limits, structural domains, and buried calderas.

With respect to (7), it is important to realize that SP surveys do not directly image these limits and structural features. SP signals highlight their impact on groundwater flow and hydrothermal circulation, from which structural inferences can be made.

Through integration with different geophysical and geochemical data types such as subsurface electrical resistivity distributions obtained from, e.g., electrical resistivity tomography or magnetotellurics, soil CO₂ flux, and soil temperature, the hydrogeological interpretations obtained from SP measurements can be better constrained and/or validated. Furthermore, knowing the subsurface electrical resistivity distribution, SP data can be used to estimate Darcy flow velocities that can directly inform hydrologic models [44].

However, groundwater and hydrothermal fluid flow is inherently a 3D problem.

Using networks of 2D ERT acquisition lines, a pseudo-3D electrical conductivity image of the subsurface can be obtained [91]. With emerging nodal-based electrical resistivity tomography and Induced Polarization systems that measure both horizontal components of the electric field, such as the pole-dipole FullWaver-V systems from Iris Instruments, we can obtain true 3D electrical resistivity and induced polarization images, benefiting from the fact that the system allows for flexible acquisition geometries, data acquisition in areas with steep topography, and offers much larger depths of investigation compared to cabled ERT-systems (due to the FullWaver-V pole-dipole configuration). Combining these 3D datasets with 2D SP maps will enable studying groundwater and hydrothermal fluid flow in 3D, which will enhance our understanding of groundwater and hydrothermal systems and will be highly valuable for groundwater modelers.

Alternatively, for deeper targets, 3D resistivity models can be obtained through, e.g., 3D MT surveys [32]. Traditionally, gravity surveys have been proven highly useful to obtain 3D density structural information using relatively fast and flexible acquisition, for example, in a tuff ring complex [92], for inference of the density structure of Piton de la Fournaise and mass transfer in this volcanic system [93], or for a volcanic field in Saudi Arabia [94].

In recent years, airborne geophysical surveys, for example, airborne electromagnetic surveys, have been proven very successful for imaging the groundwater system with large spatial coverage [34–37].

Furthermore, the latest developments in muon tomography demonstrate that muon imaging is a powerful technique, especially for complementing other geophysical methods: for example combined with gravity data to image the 3D density structure of lava domes [95], or together with seismic data to detect abrupt changes of hydrothermal activity [96].

Integrating SP data with any of these geophysical data types will enable a 3D understanding of groundwater flow patterns, hydrothermal systems, and volcanic structural features that control fluid flow.

For the case studies presented in this paper, the measured SP signals are interpreted as being the result of macroscale structural features and fluid flow patterns. As discussed earlier, SP signals have a microscale origin, and even though we have interpreted the case studies in this paper to be primarily the consequence of an electrokinetic mechanism, it is important to keep in mind that other microscale mechanisms may be at play (such as thermoelectric and electrochemical effects). This may be important especially in volcanic environments, where, for example, in hydrothermal systems, hot, acidic fluids circulate.

Moreover, in fractured, volcanic rock formations, heterogeneities that are smaller than the characteristic scale of the problem (such as the size of typical seismic wavelengths), but larger than the pore scale, can affect fluid flow at the so-called mesoscale [97]. The mesoscale forms the bridge between the macroscale features we typically interpret using SP signals and the microscale origin of the SP signals. For seismoelectric signals, where SP-like signals are actively generated by seismic wave-induced fluid flow, it has been shown that features such as mesoscopic fractures can create measurable seismoelectric converted signals under laboratory conditions [98]. Using an effective field theory, it has furthermore been shown that three effective properties of inhomogeneous media playing a role in electrokinetic phenomena (permeability, electroosmotic coupling coefficient, and electrical conductivity) are not completely independent variables when considering mesoscale inhomogeneities [99].

However, an interpretation of SP signals on the mesoscale is not often made. Given the strong macroscale imprint on the measured SP signals, it is likely that mesoscopic (and microscopic) effects are often hard to recover. However, for example, in scenarios of strong preferential fluid flow through dense networks of microcracks and mesoscale fractures, which can be highly probable in volcanic environments, it is very likely that a mesoscale formulation of the SP problem and accompanying interpretation will become important. Future research in this direction will be carried out.

Despite the fact that microscale and mesoscale properties have not been inferred in this paper, the case studies clearly show that, using the self-potential method, we have a cheap and efficient way of studying volcanic systems in terms of their macroscale structural features and flow paths, regardless of the type of volcanic environment (e.g., tuff cones, shield volcanoes, stratovolcanoes, and monogenetic fields).

Given the demonstrated efficiency, success, and value of self-potential measurements to study the hydrogeology, as well as the hydrologic and hydrothermal processes of volcanic environments, we expect SP-measurements to be highly valuable in our current, multigeophysical efforts as part of the multidisciplinary NSF EPSCoR ‘Ike Wai project, aiming to better understand the hydrology of the Hawaiian Islands for a sustainable water future.

Data Availability

The data used in this paper are available upon request to the authors.

Disclosure

S. Barde-Cabusson is currently 100% funded by the NSF EPSCoR ‘Ike Wai project, and N. Grobbe is funded for 50% by the NSF EPSCoR ‘Ike Wai project and for 50% by the University of Hawai‘i at Mānoa Water Resources Research Center. Support for the Hawai‘i EPSCoR Program is provided by the National Science Foundation’s Research Infrastructure Improvement (RII) Track-1: ‘Ike Wai: Securing Hawaii's Water Future Award # OIA-1557349.

Conflicts of Interest

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

Acknowledgments

The authors would like to explicitly acknowledge the various scientists and colleagues whose time, energy, and scientific contributions resulted in the case studies displayed in this paper. N. Grobbe would like to warmly thank Steve Martel for inspiring discussions on fractures. The authors thank the editor and an anonymous reviewer for the great suggestions and comments that helped to improve this manuscript. This is contributed paper #2394 of the Hawai‘i Institute of Geophysics and Planetology, #10765 of the School of Ocean and Earth Science and Technology, and WRRC-CP-2020-03 of the Water Resources Research Center; University of Hawai‘i at Manoa, Honolulu, Hawai’i, USA.

References

K. Takasaki, G. Hirashima, and E. Lubke, “Water resources of windward o‘ahu, hawaii,” Geological Survey Water Supply Paper, vol. 1894, pp. 1–119, 1969.

S. Lau and J. Mink, Hydrology of the Hawaiian Islands, University of Hawaii Press, Honolulu, Hawaii, USA, 2006.

N. Lautze, D. Thomas, N. Hinz, G. Apuzen-Ito, N. Frazer, and D. Waller, “Play fairway analysis of geothermal resources across the State of Hawaii: 1. Geological, geophysical, and geochemical datasets,” Geothermics, vol. 70, pp. 376–392, 2017.

A. Soueid Ahmed, A. Revil, S. Byrdina et al., “3D electrical conductivity tomography of volcanoes,” Journal of Volcanology and Geothermal Research, vol. 356, pp. 243–263, 2018.

A. Finizola, T. Ricci, R. Deiana et al., “Adventive hydrothermal circulation on Stromboli volcano (Aeolian Islands, Italy) revealed by geophysical and geochemical approaches: Implications for general fluid flow models on volcanoes,” Journal of Volcanology and Geothermal Research, vol. 196, no. 1-2, pp. 111–119, 2010.

E. Anderson, The Dynamics of Faulting and Dike Formation, Oliver and Boyd, London, UK, 1951.

H. Odé, “Mechanical analysis of the dike pattern of the Spanish peaks area, Colorado,” Bulletin of the Geological Society of America, vol. 68, no. 5, pp. 567–576, 1957.

L. G. Mastin and D. D. Pollard, “Surface deformation and shallow dike intrusion processes at Inyo craters, Long Valley, California,” Journal of Geophysical Research: Atmospheres, vol. 93, no. 11, pp. 13221–13235, 1988.

K. A. Daniels, J. L. Kavanagh, T. Menand, and R. S. J. Sparks, “The shapes of dikes: Evidence for the influence of cooling and inelastic deformation,” Bulletin of the Geological Society of America, vol. 124, no. 7-8, pp. 1102–1112, 2012.

M. Holland, J. L. Urai, and S. Martel, “The internal structure of fault zones in basaltic sequences,” Earth and Planetary Science Letters, vol. 248, no. 1-2, pp. 301–315, 2006.

G. Müller, “Experimental simulation of basalt columns,” Journal of Volcanology and Geothermal Research, vol. 86, no. 1-4, pp. 93–96, 1998.

M. P. Ryan and C. G. Sammis, “Cyclic fracture mechanisms in cooling basalt,” Bulletin of the Geological Society of America, vol. 89, no. 9, pp. 1295–1308, 1978.

C. Carrigan, “A two-phase hydrothermal cooling model for shallow intrusions,” Journal of Volcanology and Geothermal Research, vol. 28, no. 1-2, pp. 175–192, 1986.

S. Martel, “Mechanics of landslide initiation as a shear fracture phenomenon,” Marine Geology, vol. 203, no. 3-4, pp. 319–339, 2004.

R. M. Iverson, “Can magma-injection and groundwater forces cause massive landslides on Hawaiian volcanoes?” Journal of Volcanology and Geothermal Research, vol. 66, no. 1-4, pp. 295–308, 1995.

E. Custodio, “Groundwater in volcanic hard rocks,” in Groundwater in Fractured Rocks, pp. 95–108, Taylor & Francis, Milton Park, UK, 2007.

N. Van Breemen and P. Buurman, “Ferralitisation,” in Soil Formation, Chapter 13, pp. 307–328, Springer, Dordrecht, 2003.

P. Lachassagne, B. Aunay, N. Frissant, M. Guilbert, and A. Malard, “High-resolution conceptual hydrogeological model of complex basaltic volcanic islands: a Mayotte, Comoros, case study,” Terra Nova, vol. 26, no. 4, pp. 307–321, 2014.

C. Herrera and E. Custodio, “Conceptual hydrogeological model of volcanic Easter Island (Chile) after chemical and isotopic surveys,” Hydrogeology Journal, vol. 16, no. 7, pp. 1329–1348, 2008.

A. Revil and A. Jardani, The Self-Potential Method - Theory and Applications in Environmental Geosciences, Cambridge University Press, New York, NY, USA, 2013.

J. Join, J. Folio, and B. Robineau, “Aquifers and groundwater within active shield volcanoes. Evolution of conceptual models in the Piton de la Fournaise volcano,” Journal of Volcanology and Geothermal Research, vol. 147, no. 1-2, pp. 187–201, 2005.

S. Planke, E. Alvestad, and O. Eldholm, “Seismic characteristics of basaltic extrusive and intrusive rocks,” The Leading Edge, vol. 18, no. 3, pp. 342–348, 1999.

F. Brenguier, N. M. Shapiro, M. Campillo, A. Nercessian, and V. Ferrazzini, “3-D surface wave tomography of the Piton de la Fournaise volcano using seismic noise correlations,” Geophysical Research Letters, vol. 34, pp. 1–5, 2007.

A. Mordret, D. Rivet, M. Landès, and N. M. Shapiro, “Three-dimensional shear velocity anisotropic model of Piton de la Fournaise Volcano (La Réunion Island) from ambient seismic noise,” Journal of Geophysical Research: Solid Earth, vol. 120, no. 1, pp. 406–427, 2015.

E. Larose, S. Carrière, C. Voisin et al., “Environmental seismology: What can we learn on earth surface processes with ambient noise?” Journal of Applied Geophysics, vol. 116, pp. 62–74, 2015.

H. Yao and R. D. van der Hilst, “Analysis of ambient noise energy distribution and phase velocity bias in ambient noise tomography, with application to SE Tibet,” Geophysical Journal International, vol. 179, no. 2, pp. 1113–1132, 2009.

J. Robinson, L. Slater, T. Johnson et al., “Imaging pathways in fractured rock using three-dimensional electrical resistivity tomography,” Groundwater, vol. 54, no. 2, pp. 186–201, 2016.

J.-F. Lénat, D. Fitterman, D. B. Jackson, and P. Labazuy, “Geoelectrical structure of the central zone of Piton de la Fournaise volcano (Reunion),” Bulletin of Volcanology, vol. 62, no. 2, pp. 75–89, 2000.

A. Kaufman and G. Keller, Frequency and Transient Sounding, Elsevier, Amsterdam, The Netherlands, 1983.

E. M. Stolper, D. J. DePaolo, and D. M. Thomas, “Deep drilling into a mantle plume volcano: the hawaii scientific drilling project,” Scientific Drilling, no. 7, pp. 4–14, 2009.

H. Pierce and D. Thomas, “Magnetotelluric and audiomagnetotelluric groundwater survey along the Humu'ula portion of Saddle Road near and around the Pohakuloa Training Area, Hawaii,” Open-File Report 2009-1135, Hawaii, 2009.

M. Rosas-Carbajal, N. Linde, J. Peacock, F. Zyserman, T. Kalscheuer, and S. Thiel, “Probabilistic 3-D time-lapse inversion of magnetotelluric data: application to an enhanced geothermal system,” Geophysical Journal International, vol. 203, no. 3, pp. 1946–1960, 2015.

C. Hogg, D. Kiyan, V. Rath et al., “3-D interpretation of short-period magnetotelluric data at Furnas Volcano, Azores Islands,” Geophysical Journal International, vol. 213, no. 1, pp. 371–386, 2018.

A. Pryet, N. D'Ozouville, S. Violette, B. Deffontaines, and E. Auken, “Hydrogeological settings of a volcanic island (San Cristóbal, Galapagos) from joint interpretation of airborne electromagnetics and geomorphological observations,” Hydrology and Earth System Sciences, vol. 16, no. 12, pp. 4571–4579, 2012.

M. Dumont, P. Reninger, A. Pryet, G. Martelet, B. Aunay, and J. Join, “Agglomerative hierarchical clustering of airborne electromagnetic data for multi-scale geological studies,” Journal of Applied Geophysics, vol. 157, pp. 1–9, 2018.

B. Vittecoq, P. Reninger, F. Lacquement, G. Martelet, and S. Violette, “Hydrogeological conceptual model of andesitic watersheds revealed by high-resolution heliborne geophysics,” Hydrology and Earth System Sciences, vol. 23, no. 5, pp. 2321–2338, 2019.

R. Knight, R. Smith, T. Asch et al., “Mapping aquifer systems with airborne electromagnetics in the central valley of california,” Groundwater, vol. 56, no. 6, pp. 893–908, 2018.

D. B. Jackson and J. Kauahikaua, “Regional self-potential anomalies at Kilauea Volcano (Hawaii),” US Geological Survey Professional Paper, vol. 1350, no. 2, pp. 947–959, 1987.

T. Hashimoto and Y. Tanaka, “A large self-potential anomaly on Unzen Volcano, Shimabara Peninsula, Kyushu Island, Japan,” Geophysical Research Letters, vol. 22, no. 3, pp. 191–194, 1995.

A. Finizola, F. Sortino, J. Lénat, and M. Valenza, “Fluid circulation at Stromboli volcano (Aeolian Islands, Italy) from self-potential and CO2 surveys,” Journal of Volcanology and Geothermal Research, vol. 116, no. 1-2, pp. 1–18, 2002.

T. Ishido, “Electrokinetic mechanism for the “W”-shaped self-potential profile on volcanoes,” Geophysical Research Letters, vol. 31, pp. 1–5, 2004.

A. Finizola, A. Revil, E. Rizzo et al., “Hydrogeological insights at Stromboli volcano (Italy) from geoelectrical, temperature, and CO2 soil degassing investigations,” Geophysical Research Letters, vol. 33, no. 17, pp. 1–4, 2006.

A. Bolève, A. Revil, F. Janod, J. L. Mattiuzzo, and J.-J. Fry, “Preferential fluid flow pathways in embankment dams imaged by self-potential tomography,” Near Surface Geophysics, vol. 7, no. 5-6, pp. 447–462, 2009.

A. Jardani, A. Revil, A. Bolève et al., “Tomography of the Darcy velocity from self-potential measurements,” Geophysical Research Letters, vol. 34, no. 24, pp. 1–6, 2007.

G. Petiau and A. Dupis, “Noise, temperature coefficient and long time stability of electrodes for telluric observations,” Geophysical Prospecting, vol. 28, no. 5, pp. 792–804, 1980.

G. Petiau, “Second generation of lead-lead chloride electrodes for geophysical applications,” Pure and Applied Geophysics, vol. 157, no. 3, pp. 357–382, 2000.

J. E. Nyquist and C. E. Corry, “Self-potential: The ugly duckling of environmental geophysics,” The Leading Edge, vol. 21, no. 5, pp. 446–451, 2002.

F. Schoemaker, Electrokinetic Conversion, Delft University of Technology: Proefschriftmaken.nl, The Netherlands, 2011.

L. Onsager, “Reciprocal relations in irreversible processes- I,” Physical Review A: Atomic, Molecular and Optical Physics, vol. 37, no. 4, pp. 405–426, 1931.

L. Onsager, “Reciprocal relations in irreversible processes- II,” Physical Review A: Atomic, Molecular and Optical Physics, vol. 38, no. 12, pp. 2265–2279, 1931.

O. Stern, “The theory of the electrolytic double-layer,” Elektro-chemistry, vol. 30, article 508, 1924.

M. Gouy, “Sur la constitution de la électrique a la surface d'un electrolyte,” Journal de Physique Théorique et Appliquée, vol. 9, no. 1, pp. 457–468, 1910.

D. Chapman, “A contribution to the theory of electrocapillarity,” Philosophical Magazine, vol. 25, pp. 475–481, 1913.

A. Revil, A. Jardani, P. Sava, and A. Haas, The Seismoelectric Method - Theory and Applications, John Wiley & Sons, Chichester, England, 2015.

B. Minsley, Modeling and Inversion of Self-Potential Data, MIT, Cambridge, UK, 2007.

L. Jouniaux, A. Maineult, V. Naudet, M. Pessel, and P. Sailhac, “Review of self-potential methods in hydrogeophysics,” Comptes Rendus Geoscience, vol. 341, no. 10-11, pp. 928–936, 2009.

E. Leinov, J. Vinogradov, and M. D. Jackson, “Salinity dependence of the thermoelectric coupling coefficient in brine-saturated sandstones,” Geophysical Research Letters, vol. 37, pp. 1–6, 2010.

A. Revil, C. D. Meyer, and Q. Niu, “A laboratory investigation of the thermoelectric effect,” Geophysics, vol. 81, no. 4, pp. E243–E257, 2016.

A. Revil, “Ionic diffusivity, electrical conductivity, membrane and thermoelectric potentials in colloids and granular porous media: a unified model,” Journal of Colloid and Interface Science, vol. 212, no. 2, pp. 503–522, 1999.

L. Jouniaux and T. Ishido, “Electrokinetics in earth sciences: a tutorial,” International Journal of Geophysics, vol. 2012, Article ID 286107, 16 pages, 2012.

J. Zlotnicki and Y. Nishida, “Review on morphological insights of self-potential anomalies on volcanoes,” Surveys in Geophysics, vol. 24, no. 4, pp. 291–338, 2003.

M. Aubert, I. Dana, and M. Livet, “Vérification de limites de nappes aquifères en terrain volcanique par la méthode de polarisation spontanée,” Les Comptes Rendus de l'Académie des Sciences, vol. 311, no. II, pp. 999–1004, 1990.

S. Barde-Cabusson, A. Finizola, A. Revil et al., “New geological insights and structural control on fluid circulation in La Fossa cone (Vulcano, Aeolian Islands, Italy),” Journal of Volcanology and Geothermal Research, vol. 185, no. 3, pp. 231–245, 2009.

A. Revil, A. Finizola, S. Piscitelli et al., “Inner structure of La Fossa di Vulcano (Vulcano Island, southern Tyrrhenian Sea, Italy) revealed by high-resolution electric resistivity tomography coupled with self-potential, temperature, and CO2 diffuse degassing measurements,” Journal of Geophysical Research: Solid Earth, vol. 113, pp. 1–21, 2008.

G. Frazzetta, L. La Volpe, and M. F. Sheridan, “Evolution of the Fossa Cone, Vulcano,” Journal of Volcanology and Geothermal Research, vol. 17, no. 1-4, pp. 329–360, 1983.

P. Dellino and L. La Volpe, “Stratigrafia, dinamiche eruttive e deposizionali, scenario eruttivo e valutazioni di pericolositα a La Fossa di vulcano,” La Volpe, Felici, Pisa., CNR-GNV Progetto Vulcano 1993-1995, pp. 214–237, 1997.

R. De Rosa, N. Calanchi, P. Dellino et al., “Field Trip Guide Book - Geology and volcanism of Stromboli, Lipari, and Vulcano,” in Proceedings of the in 32nd International Geological Congress, Firenze, 2004.

G. De Astis, P. Dellino, L. La Volpe, F. Lucchi, and C. Tranne, Geological map of the Vulcano Island 1:10000, 2007.

J.-F. Lénat and P. Bachèlery, “Structure et fonctionnement de la zone centrale du Piton de la Fournaisse,” in Le volcanisme de la Réunion - Monographie, Clermond-Ferrand, France, Cent. Rech. Volcanol, pp. 257–296, 1990.

L. Michon, T. Staudacher, V. Ferrazzini, P. Bachèlery, and J. Marti, “April 2007 collapse of Piton de la Fournaise: A new example of caldera formation,” Geophysical Research Letters, vol. 34, no. L21301, pp. 1–4, 2007.

M. Urai, N. Geshi, and T. Staudacher, “Size and volume evaluation of the caldera collapse on Piton de la Fournaise volcano during the April 2007 eruption using ASTER stereo imagery,” Geophysical Research Letters, vol. L22318, pp. 1–4, 2007.

S. Barde-Cabusson, A. Finizola, A. Peltier et al., “Structural control of collapse events inferred by self-potential mapping on the Piton de la Fournaise volcano (La Réunion Island),” Journal of Volcanology and Geothermal Research, vol. 209, pp. 9–18, 2012.

A. Abdel-Monem, N. D. Watkins, and P. W. Gast, “Potassium-argon ages, volcanic stratigraphy, and geomagnetic polarity history of the Canary Islands; Tenerife, La Palma and Hierro,” American Journal of Science, vol. 272, no. 9, pp. 805–825, 1972.

E. Ancochea, J. Fuster, E. Ibarrola et al., “Volcanic evolution of the island of Tenerife (Canary Islands) in the light of new K-Ar data,” Journal of Volcanology and Geothermal Research, vol. 44, no. 3-4, pp. 231–249, 1990.

J. Marti and A. Gudmundsson, “The Las Cañadas caldera (Tenerife, Canary Islands): an overlapping collapse caldera generated by magma-chamber migration,” Journal of Volcanology and Geothermal Research, vol. 103, no. 1-4, pp. 161–173, 2000.

V. Villasante-Marcos, A. Finizola, R. Abella et al., “Hydrothermal system of Central Tenerife Volcanic Complex, Canary Islands (Spain), inferred from self-potential measurements,” Journal of Volcanology and Geothermal Research, vol. 272, pp. 59–77, 2014.

J. Martí, A. Geyer, J. Andujar, F. Teixidó, and F. Costa, “Assessing the potential for future explosive activity from Teide–Pico Viejo stratovolcanoes (Tenerife, Canary Islands),” Journal of Volcanology and Geothermal Research, vol. 178, no. 3, pp. 529–542, 2008.

J. C. Carracedo, E. Rodríguez Badiola, H. Guillou et al., “Eruptive and structural history of Teide Volcano and rift zones of Tenerife, Canary Islands,” Bulletin of the Geological Society of America, vol. 119, no. 9-10, pp. 1027–1051, 2007.

I. Domínguez-Cerdeña, C. del Fresno, and L. Rivera, “New insight on the increasing seismicity during Tenerife's 2004 volcanic reactivation,” Journal of Volcanology and Geothermal Research, vol. 206, no. 1-2, pp. 15–29, 2011.

S. Barde-Cabusson, G. Levieux, J. Lénat et al., “Transient self-potential anomalies associated with recent lava flows at Piton de la Fournaise volcano (Réunion Island, Indian Ocean),” Journal of Volcanology and Geothermal Research, vol. 187, no. 3-4, pp. 158–166, 2009.

M. Rosi, “The island of Stromboli,” Rendiconti della Societa Italiana di Mineralogia e Petrologia, vol. 36, pp. 345–368, 1980.

G. Pasquarè, L. Francalanci, V. Garduno, and A. Tibaldi, “Structure and geologic evolution of the Stromboli volcano, Aeolian Islands, Italy,” Acta Vulcanologica, vol. 3, pp. 79–89, 1993.

F. Barberi, M. Rosi, and A. Sodi, “Volcanic hazard assessment at Stromboli based on review of historical data,” Acta Vulcanol, vol. 3, pp. 173–187, 1993.

N. Coppo, P. Schnegg, W. Heise, P. Falco, and R. Costa, “Multiple caldera collapses inferred from the shallow electrical resistivity signature of the Las Cañadas caldera, Tenerife, Canary Islands,” Journal of Volcanology and Geothermal Research, vol. 170, no. 3-4, pp. 153–166, 2008.

G. Guérin and G. Valladas, “Thermoluminescence dating of volcanic plagioclases,” Nature, vol. 286, no. 5774, pp. 697–699, 1980.

G. Guérin, G. Behamoun, J. Mallarach, and G. Guérin, “Un exemple de fusió parcial en medi continental. El vulcanisme quaternari de la Garrotxa,” Publicació del Museu Comarcal de la Garrotxa, vol. 1, pp. 19–26, 1985.

J. Cebriá, J. López-Ruiz, M. Doblas, R. Oyarzun, J. Hertogen, and R. Benito, “Geochemistry of the Quaternary alkali basalts of Garrotxa (NE Volcanic Province, Spain): a case of double enrichment of the mantle lithosphere,” Journal of Volcanology and Geothermal Research, vol. 112, no. 3-4, pp. 175–187, 2000.

X. Bolós, S. Barde-Cabusson, D. Pedrazzi et al., “Geophysical exploration on the subsurface geology of La Garrotxa monogenetic volcanic field (NE Iberian Peninsula),” International Journal of Earth Sciences, pp. 1–15, 2014.

S. Barde-Cabusson, J. Gottsmann, J. Martí et al., “Structural control of monogenetic volcanism in the Garrotxa volcanic field (Northeastern Spain) from gravity and self-potential measurements,” Bulletin of Volcanology, vol. 76, no. 788, pp. 1–13, 2014.

J. Lénat, “Retrieving self-potential anomalies in a complex volcanic environment: an SP/elevation gradient approach,” Near Surface Geophysics, vol. 5, no. 3, pp. 161–170, 2007.

M. Rosas-Carbajal, J.-C. Komorowski, F. Nicollin, and D. Gibert, “Volcano electrical tomography unveils edifice collapse hazard linked to hydrothermal system structure and dynamics,” Scientific Reports, vol. 6, pp. 1–11, 2016.

G. Ito, J. Higa, R. Whittier, N. Grobbe, N. Lautze, and D. Thomas, “Constraints on subsurface density from gravity surveying in and around the salt lake tuff ring complex, south-central oahu,” SOEST Technical Report 10644, Honolulu, HI, USA, 2019.

L. Gailler, J. Lénat, M. Lambert, G. Levieux, N. Villeneuve, and J. Froger, “Gravity structure of Piton de la Fournaise volcano and inferred mass transfer during the 2007 crisis,” Journal of Volcanology and Geothermal Research, vol. 184, no. 1-2, pp. 31–48, 2009.

E. Aboud, A. Shareef, F. A. Alqahtani, and S. Mogren, “Using a 3D gravity inversion technique to image the subsurface density structure in the Lunayyir volcanic field, Saudi Arabia,” Journal of Asian Earth Sciences, vol. 161, pp. 14–24, 2018.

M. Rosas-Carbajal, K. Jourde, J. Marteau, S. Deroussi, J.-C. Komorowski, and D. Gibert, “Three-dimensional density structure of La Soufrière de Guadeloupe lava dome from simultaneous muon radiographies and gravity data,” Geophysical Research Letters, vol. 44, no. 13, pp. 6743–6751, 2017.

Y. Le Gonidec, M. Rosas-Carbajal, J. d. Bremond d’Ars et al., “Abrupt changes of hydrothermal activity in a lava dome detected by combined seismic and muon monitoring,” Scientific Reports, vol. 9, pp. 1–9, 2019.

T. M. Müller, B. Gurevich, and M. Lebedev, “Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks - A review,” Geophysics, vol. 75, no. 5, pp. 75A147–75A164, 2010.

D. Jougnot, J. G. Rubino, M. R. Carbajal, N. Linde, and K. Holliger, “Seismoelectric effects due to mesoscopic heterogeneities,” Geophysical Research Letters, vol. 40, no. 10, pp. 2033–2037, 2013.

V. Levin, M. Markov, and I. Markova, “Effective electrokinetic parameters of inhomogeneous porous media: the effective field approximation,” Transport in Porous Media, vol. 121, no. 2, pp. 459–477, 2018.