Assessing hydrofacies and hydraulic properties of basaltic aquifers derived from geophysical logging

Basaltic aquifers are an important source of water supply in many regions worldwide. Because of their cooling process, flood basalts normally have complex internal structures and unpredictable permeable zone distribution. Geophysical profiling is a reliable tool used to identify the vertical variations of physical properties within the flood basalts, mainly to distinguish between low and high permeability intervals. In this study, a detailed analysis of drill cuttings, and geophysical logging of two depth wells helped to identify the typical response of main facies of flood basalts in the Serra Geral aquifer system. The evidence obtained from this work confirms that flood basalts may be classified as multilayer aquifer systems in which a highly permeable top covers a center with low permeability as a result of recurrent floods. In several cases, permeable layers are represented by a weathering horizon developed between two flood events. This study uses an empirical model to estimate the porosity based on acoustic velocity logging. A useful model based on a well-established Kozeny-Carman model was developed to predict the permeability of basaltic rocks using porosity data. The results obtained allow to identify the permeable intervals of basaltic aquifers and to estimate their hydraulic properties.


INTRODUCTION
Basaltic aquifers represent an important source of water, mainly for agricultural purposes, in many regions worldwide such as western India (Limaye 2010), the northwestern United States (Piersol and Sprenke 2015), Mexico City (Edmunds et al. 2002), northeastern Australia (Locsey and Cox 2003), and southern Ethiopia (McKenzie et al. 2001). Moreover, groundwater obtained from basalts represents the main source of water supply in volcanic islands such as Jeju Island, Korea (El-Kadi et al. 2014), the Canary Islands, Spain (Cabrera and Custodio 2004), Madeira Island, Portugal (Prada et al. 2005), and Hawaii, the United States (Wentworth 1951). Despite their importance, the distribution of permeable zones within flood basalts is poorly understood due to the internal complexity of these aquifers.
The distribution of permeable zones within basalts is intrinsically related to the internal structure of basalts that were generated by the cooling processes after a magma flood. The cooling structures have been exhaustively investigated in several studies such as those by James (1920), Tomkeieff (1940), Walker (1971), Long and Wood (1986), Lore et al. (2000), Kattenhorn and Schaefer (2008), Browning et al. (2016), and Puffer et al. (2018). The cooling structures create internal layers in response to distinct cooling rates within the flood, as described in Uhl and Joshi (1986), Whitehead (1992), Kulkarni et al. (2000), and McGrail et al. (2006). According to these studies, basalts typically display the same facies succession, consisting of three distinct units: flow top, flow interior, and flow bottom.
Flow tops (hereafter referred to as porous flow tops) are usually characterized by a genetic relationship between degassing, thermal contractions, and interactions with water. The tops of flood basalts are comprised of a chilled, glassy crust that is often vesicular to scoriaceous or brecciated rubble produced by escaping gas during the cooling process (Buckley and Oliver 1990, Whitehead 1992, McGrail et al. 2006, Zakharova et al. 2012. The vesicles represent the primary porosity of the rock and are generated from gas flows during lava flooding. The porosity is dependent on bubble gas and lava flow characteristics (Sahagian et al. 1989, Mangan and Cashman 1996, Melnik 2000 and can vary widely, as demonstrated through several experimental studies. Another characteristic identified by Buckley and Oliver (1990) described the flow top in the Deccan traps as soft and weathered, consisting of a horizon composed of green, grey, or red clay with amygdaloidal material, with a typical thickness of 1.5-3.0 m. McGrail et al. (2006) presented evidence of a high-lateral continuity of vesicular layers that serve as regional aquifer systems. The slow cooling within the flooded lava allows massive basalts to crystalize and become very thick. During the cooling process, contraction stress breaks the rock, generating typical colonnade and entablature structures. Long and Wood (1986) describe the colonnade as columnar structures oriented perpendicularly to flow boundaries, and typically occurring in the basal portion of the flow, alternating with entablatures. Similar to flow top, the flow bottom is comprised of a glassy, vesicular zone (e.g., Whitehead 1992, Uhl and Joshi 1986, Buckley and Oliver 1990, Zakharova et al. 2012. However, as noted by Buckley and Oliver (1990), this zone is typically only a few centimeters thick.
The facies that comprise basalt floods present strong variations of physical properties (e.g., Peterson and Lao 1970, Helm-Clark et al. 2004, Asfahani et al. 2011, Zakharova et al. 2012. The strong contrast in physical properties means that geophysical logging is a prominent tool that can be used to characterize the stratigraphic features of flood basalts (Helm-Clark et al. 2004). An extensive revision of geophysical logging was described in detail by Helm-Clark et al. (2004), while Zakharova et al. (2012) analyzed the hydraulics of core samples and geophysical logging data from basalt formations of the Columbia River Basalt Group. The porosity, represented by vesicles and fractures, is one of the main features that govern the variation of physical properties of basalts, as demonstrated by several authors such as Helm-Clark et al. (2004), Asfahani (2011), andZakharova et al. (2012).
Despite their global abundance and importance, research on basaltic aquifers is less common when compared to typical granular aquifers. Since basalt acts as an important reservoir, this study aims to identify the parameters derived from geophysical logging that allow us to diagnose the facies variation within flood basalts. Furthermore, this research aims to develop empirical models that produce reliable estimates of porosity and permeability of rocks using geophysical logging.

Study area
The Serra Geral formation represents one of the largest continental flood basalt occurrences, covering nearly 917,000 km 2 (Frank et al. 2009). According to Almeida (1986), the Serra Geral formation represents a thick sequence flow of mainly basaltic rocks, with a maximum thickness of 1,700 m, belonging to the continental-scale Paraná-Etendeka Magmatic Province. According to Turner et al. (1994), the age of magmatism ranges from 137 to 127 Ma. The stratigraphy, geochemistry, and tectonic history of Serra Geral were exhaustively studied in works including Frank et al. (2009), Machado et al. (2015, Barreto et al. (2016), Polo et al. (2018), Rossetti et al. (2018), and Svensen et al. (2017).
Wireline logs for the open section of two boreholes (340-1,250 m) were acquired in the county of Colômbia, São Paulo State, located in the southeast region of Brazil (Fig. 1). The analyzed basalt was covered by sandstone and claystone, belonging to the Bauru Group of the Cretaceous age, that were nearly 100 meters thick. The basalt also covers the sandstone of the Pirambóia and Botucatu formations.

Geophysical profile analysis
Based on descriptions provided in previous works, flow tops were interpreted as permeable zones, characterized by low resistivity, low sonic velocity, and high gamma rays. Using this log analysis coupled with rock samples, the typical basalt succession found is comprised of low permeability levels being vertically succeeded by high permeability horizons. The large intervals of geophysical profiling provide the possibility of identifying the limits of flood basalts and the proportion of permeable layers within the individual flood basalt flows. The wireline logs were conducted prior to well completion. Caliper (borehole diameter), total natural gamma radioactivity, electrical resistivity, and sonic velocities data were acquired using Schlumberger profiling tools. The drill cuttings description was also employed during well logs interpretation. Our analysis was restricted to total natural gamma radioactivity, electrical resistivity, and sonic velocities, since these combined parameters sufficiently describe the variation of flood basalts' internal structures.
Resistivity and sonic logs represent the most useful methods for identifying flow tops and flow interiors because of the strong variation in porosity observed therein. According to Helm-Clark et al. (2004), resistivity represents the best logging method to discriminate between permeable and impermeable zones of basaltic aquifers. Basalt is one of the most resistive rocks in the Earth's crust (Helm-Clark et al. 2004) and the resistivity of low-porosity flow interior is markedly high. Thus, an increase of porosity and, consequently, stored water content, promote a sharp decrease in resistivity. Porosity neutrons and sonic logs are very suitable for the identification of permeable and impermeable zone limits (Buckley and Oliver 1990, Helm-Clark et al. 2004, Zakharova et al. 2012, since porosity and permeability are correlated. The sonic P wave velocity log (Vp) corresponds to the propagation velocity of an elastic wave through the tested geological unit. Due to the differences in elastic properties in response to porosity (Zamora et al. 1994, Planke et al. 1999, Chen et al. 2015, Vp represents an important tool to estimate basalt porosity. Vp also allows to identify fractured zones (Schlumberger 1989) as indicated by an increase in the travel time required for the pulse to reach the detectors. Similar to the resistivity, acoustic velocity is high in the massive interior of flows, and low at flow tops and interbeds.
Natural rocks emit gamma rays in response to decay chains of radioisotopes 40 K, 238 U, and 232 Th and their derivative products. Thus, the gamma-ray log indirectly measures the amount of radioactive species in rocks. The gamma-ray log is normally reported in American Petroleum Institute (API) units (Schlumberger 1989). Thus, gamma-ray log is a useful tool to identify potassium-rich intervals in basalt or clay content in interbed sediments, as reported by Helm-Clark et al. (2004). Buckley and Oliver (1990) demonstrated that gamma-rays are also useful to identify weathered zones, where the abundance of clay minerals is noticeable.
Porosity and permeability derived from empirical models Sonic logs are useful techniques to estimate porosity using either the Wyllie time-average equation (Wyllie et al. 1956) or the model proposed by Raymer et al. (1980). However, these equations were developed for sedimentary rocks (mainly limestone and sandstone) and therefore have low representativeness in the estimation of basalt porosity. For a more realistic approach to quantify the porosity of basalt and to improve estimations thereof, we developed an empirical model based on sonic velocity data. To this end, experimental data relating porosity and sonic velocities in basaltic rocks, compiled from previous works (Al-Harthi et al. 1999;Zamora et al. 1994, Planke et al. 1999, Chen et al. 2015, Vedanti et al. 2018, Rossetti et al. 2019, were employed in a regression analysis. As noted by Helm-Clark et al. (2004), the use of neutron logs can potentially result in porosity overestimation in the presence of hydrous minerals, making velocity logs a better approach for porosity estimation.
The Kozeny-Carman model represents a common approach for predicting the permeability of volcanic rocks using porosity data (e.g., Klug and Cashman 1996, Saar and Manga 1999, Mueller et al. 2005, Takeuchi et al. 2008, Yokoyama and Takeuchi 2009, Degruyter et al. 2010, Farquharson et al. 2015, Lamur et al. 2017. To identify a model capable of predicting permeability from porosity, porosity and permeability measurements were obtained from previous literature. Despite the abundance of experimental data related to the porosity and permeability of general volcanic rocks, our analysis was restricted to basaltic rocks, gathering 206 pieces of experimental data relating permeability and porosity, presented by Saar and Manga (1999), Mueller et al. (2005), Loaiza et al. (2012), Lamur et al. (2017), Rossetti et al. (2019), Schaefer et al. (2015), Bai et al. (2010), and Luhmann et al. (2017). By means of non-linear regression, we built an empirical equation to estimate permeability derived from porosity values, based on the Kozeny-Carman model.

Internal structures of flood basalt
Basalts may be classified as multilayered aquifer systems, on which low-permeable levels (flow interiors) are recurrently overlapped by high-permeability zones (flow tops). This alternation of intervals is highly comparable to the genetic nature of successive covering of basaltic lava flow events. The framework of a basalt flow generated from cooling creates a vesicular and permeable zone on the flow top and a massive, entablature, or colonnade zones in the flow interior. This alternating pattern was identified in the wireline logs presented in this work and in previous ones, such as those by Buckley and Oliver (1990),  The description of drill cutting samples collected at 2 m intervals allowed us to identify 4 major facies within the drilled wells: flow interior (massive, colonnade, entablature basalts), flow tops (vesicular basalt), weathered basalt, and clay lenses. Because each facies that comprises flood basalts has a distinct range of variation in terms of density, permeability, and porosity, the responses of each geophysical method are dissimilar. Thus, the association of range values of geophysical parameters offers an unambiguous diagnostic of flood basalt internal facies. Figure 4 illustrates the boxplot of distribution of gamma-ray, electrical resistivity, and sonic P wave velocity values measured in the two wireline logs.

Flow interior
The flow interior -composed of massive basalt and/ or layers of colonnade and entablature generated by slow cooling within the center of lava flow -is characterized by low porosity and permeability and, as a result, by high resistivity (> 1,000 Ohm.m) and high Vp (mostly > 5 km/s), as seen in Figure 4. These features are unambiguously associated with massive dark gray or black basalts. A detailed description of SAG/P2 and SAG/P3 drill cuttings confirms that the electrical resistivity is the most important parameter to distinguish between high-permeability flow tops and the low-permeability flow interiors, as described by Helm-Clark et al. (2004).
Due to the remarkable contrast in resistivity, Helm-Clark et al. (2004) argue that resistivity logs represent the most useful tool to identify basalt facies, which are associated with low resistivity in interflow zones and in interbeds, and high resistivity in flow interior. Figure 4 reveals that the median resistivity of massive basalt is 563.24 Ω⋅m. Moreover, the high variability of electrical resistivity in the evaluated log, ranging from 20.6 to 100,000 Ω⋅m represents a wider range than that recorded by Zakharova et al. (2012) (i.e., 10 to 2,000 Ω⋅m). Helm-Clark et al. (2004) describe that the increase in acoustic velocity is associated with massive basalt, while the decrease of acoustic velocity may be associated with sedimentary interbeds, fractured zones, or vesicular basalts. Higher Vp is directly related to both higher density and low-porosity of massive basalt. According to the experimental verification by Zamora et al. (1994), Al-Harthi et al. (1999), Planke et al. (1999) and others, the travel time of acoustic P wave velocity is intrinsically related to basalt porosity. The median value of Vp associated with flow interior is 5.6 km/s (Fig. 4), noticeably above the values of vesicular, weathered, and clay horizons. The range of sonic velocities observed from massive basalts in the present study (4.03-5.98 km/s) is similar to that presented by Zakharova et al. (2012) in the Columbia River Basalt Group (CRBG).

Vesicular basalt
Contrary to flow interior, the vesicular basalt is notably marked by its rapid cooling. The escape of bubble gas during cooling creates a network of interconnected pores. As illustrated in Figures 2 and 3, the vesicular zones comprise a limited proportion of a basalt flow. The most noticeable difference was observed in the electrical resistivity of vesicular intervals. While Zakharova et al. (2012) recorded a narrow variation of electrical resistivity, ranging from 10-20 Ωm, our work found variations between 2.26 and 6,462.78 Ωm, with a median value of 39.27 Ωm. These dissimilarities are probably related to porosity and/or pore fluid characteristics.
The median Vp we identified was 3.9 km/s, significantly below the sonic velocity for flow interior and weathered basalts, but close to the clay lenses value. Zakharova et al. (2012) noted that the sonic velocities in vesicular basalt are about 50% of the values found in massive basalt. We found that the median Vp of vesicular and interior flow diverges around 30%. This difference may be attributed to variations of porosity between the massive and vesicular basalts of Serra Geral and Columbia River Basalts.
The vesicular intervals of basalt have high average gamma-ray values of 15.58 API, close to the values in weathered basalt and slightly higher than interior flow. Thus, gamma-rays are not useful to diagnose vesicular basalts. However, amygdaloidal basalts are strongly subjected to weathering due to the reactivity of glassy material and high-porosity, producing clays. As noted by Buckley and Oliver (1990) in Deccan traps, vesicular basalts may present high gamma-rays whenever there is pronounced weathering.

Weathered basalt
The analysis of drill cutting samples indicated that some intervals possess centimetric rounded and brown fragments, containing vesicles filled with zeolites, carbonates, and chlorite, associated with clay intervals. The development of physical and chemical weathering is likely to be intensive in basalts due to the high reactivity of major minerals in these rocks (olivines, pyroxenes, and plagioclases). Weathering progress increases the aperture of fractures and leachates the silicates, promoting a strong increase in porosity and permeability. The time interval between two lava flow events can play a crucial role in the development of an aquifer zone within the basaltic units. After lava consolidation, chemical and physical weathering take place, creating a weathering profile where permeability and porosity values are especially higher. Both Buckley and Oliver (1990) and Gingerich and Izuka (1997) verifiedon the Deccan Plateau, India and in Hawaii, the United States, respectively -the existence of weathered basalt in deep portions overlapped by non-weathered basalts. On the other hand, the core description presented by Zakharova et al. (2012) did not find weathered basalts, indicating that the weathered facies are not widespread in all basalts.
Resistivity is commonly considered the main tool to distinguish low-porosity flow interior and high-porosity flow top, since resistivity is very sensitive to changes in porosity, water-filling pores, and the presence of clay. In response to high-porosity and to the presence of clay minerals, the weathered basalt has an electrical resistivity median of 66 Ωm (Fig. 4B), a value between the flow interiors and the vesicular or clay intervals.
As demonstrated by Navarre-Sitchler et al. (2015), during the chemical weathering of basalt, rock porosity strongly increases. Thus, chemical weathering enables the use of acoustic velocity as an important tool to identify weathered intervals, since the porosity enhancement is diagnosed by the decrease in P wave acoustic velocity. The median Vp of weathered basalt is 4.5 km/s, consistently below the flow interior value (5.6 km/s), and above the sonic velocities of clay and vesicular intervals (4.3 and 3.9 km/s, respectively), as seen in Figure 4C.
Weathered basalts present a mean gamma radiation of 25.64 API, highly divergent from flow interior and vesicular basalts. The increase in gamma radiation may be related to the presence of clay minerals produced by pyroxene and plagioclase weathering. Thus, gamma radiation represents another parameter used to identify weathered basalts due to the abundance of clay minerals.

Clay lenses
Weathered horizons are associated with the presence of thick lenses of clay reported in this study as well as in previous works (e.g., Buckley and Oliver 1990, Gingerich and Izuka 1997, Jalludin and Razack 2004, and are evidence of ancient weathering surface development. Interestingly, clay lenses were not recorded in the Columbia River Basalts (Zakharova et al. 2012).
As stated by Helm-Clark et al. (2004), gamma-ray logs are unfeasible tools to identify anomalous intervals containing high potassium concentrations, which is often related to the presence of clay minerals (mainly smectite) produced by basalt weathering processes. Buckley and Oliver (1990) and Crosby and Anderson (1971) demonstrated the applicability of gamma-ray profiling tools for identification of clay lenses produced by basalt weathering. In the present study, the presence of weathered basalt associated with clay lenses was also easily identified by gamma-ray logs interpretation.
The median resistivity obtained for the clay lenses was 48.35 Ωm, close to vesicular (39.27 Ωm) and weathered (66 Ωm) intervals. The median Vp of clay lenses is 4.3 km/s, below that of vesicular basalts, and very different from weathered basalts and flow interior especially. This signature is compatible with the high-porosity of clay lenses.

Porosity and permeability derived from wireline logs
In order to build an empirical model able to produce reliable estimates of porosity from acoustic velocity, we compiled experimental data from Zamora et al. (1994), Al-Harthi et al. (1999), Planke et al. (1999), Chen et al. (2015), Vedanti et al. (2018), and Rossetti et al. (2019) relating porosity of basalts to sonic velocity (Fig. 5). Using non-linear regression, an exponential equation was fitted to 197 data points of porosity and V p , including samples spanning from massive to vesicular basalts. Due to sample size, the identified equation from non-linear regression (Eq. 1) is representative of all facies within the basalt flow. . e The porosity derived from sonic logs represents a quantitative parameter used to determine the permeable and impermeable intervals of basaltic aquifers. The contact between two successive basalt flows is characterized by an abrupt change in porosity, as seen in Figure 6. The recurrent overlap of high-(vesicular, amygdaloidal, and scoria) and low-porosity zones (massive, columnar, and entablature basalts) is evident. We estimated porosity based on Equation 1, building porosity  profiles for wells SAG/P2 and SAG/P3. In order to compare with wireline logs from the Columbia Plateau, we also derived porosity using Vp logs presented by Zakharova et al. (2012). The porosity profiles of SAG/P2, SAG/P3, and the well from Zharakova et al. (2012) are presented in Figure 6.
The porosity profile estimated by Equation 1 (Fig. 6) produces well-delimited intervals of porosity in complete agreement with the limits of basalt hydrofacies. Furthermore, Figure 6 reveals that basalts of the Columbia Plateau and the Serra Geral Aquifer have the same internal structure, displaying alternations of high and low porosity zones as well as a comparable range of porosity values. The main differences are related to the thickness of high porosity zones. While the intervals containing high-porosity in the Columbia Plateau have an average thickness of approximately 8 m, those of Serra Geral Aquifer have a thickness of approximately 12 m, reflecting the volume of lava expelled as well as its temperature and magma composition. The porosity of flow tops (vesicular, breccia, and scoria) represents the primary porosity and is generated by degassing during the cooling of sub-aerial flows, and is dependent on bubble gas and lava flow characteristics (Sahagian et al. 1989, Mangan and Cashman 1996, Melnik 2000. Vesicular basalt is characterized by high porosity due to its genetic context. The porosity of vesicular basalt was estimated using Equation 1 as ranging from 20 to 60%, agreeing with experimental results presented by Saar and Manga (1999), Song et al. (2001), Shea et al. (2010), andSchaefer et al. (2015). Reported anomalously high porosity values are associated with scoriaceous rocks, with porosity above 40% Manga 1999, Song et al. 2001), which were not found on the studied wells.
In general, the porosity of the flow interior is very low (e.g., Christensen and Wilkens 1982, Vinciguerra et al. 2005, Zakharova et al. 2012. While Vinciguerra et al. (2005) reported a porosity of 1% for columnar basalt, Christensen and Wilkens (1982) reported values in the range of 0.1-1.9%. Zakharova et al. (2012) verified that flow interiors are characterized by porosity values ranging from 0 to 10%, where higher porosity values (5-10%) correspond to higher planar fracture density. The highest porosity values of flow interiors were reported by Deolankar (1980); 15% were related to fractured basalts. Figure 6 reveals that our model predicted porosities around 4% for massive basalts in the Serra Geral Formation and the Columbia Plateau.
Experimental porosity and permeability data encompass all lithologies found in flood basalts, with porosity ranging from 1.07 to 87% and permeability ranging from 4.34 x 10 -18 to 6.89 x 10 -11 m 2 . Using non-linear regression, we established a model able to predict the permeability based on porosity data (Eq. 2). The goodness of fit was obtained using the leastsquares method of optimizing the Kozeny-Carman model by minimization of the average error (Fig. 7).
The scatterplot graph of natural logarithm of predicted permeability by Equation 2 versus observed permeability (Fig. 7B) indicates a good agreement of both values, supporting the ability of our model to produce reasonable permeability estimation. The minimum porosity value calculated in basalt was 0%, while the maximum was 39.4%. A rough estimation of permeability can be performed by applying Equation 2 and sonic-derived porosity data (Figs. 8 and 9).
In this study, the ranges of permeability values estimated (Figs. 8 and 9) produce values compatible with those recorded   in the literature. The reported permeability ranges determined for the flow tops (vesicular, breccia, and scoria) are usually high, with values ranging from 1.7 x 10 -15 to 2.4 x 10 -12 m 2 (Saar and Manga 1999, McGrail et al. 2001, Schaefer et al. 2015, Lamur et al. 2017. The permeability and porosity of the interior part of basalt flows (massive, colonnade, and entablature) are less commonly studied than those of flow tops. McGrail et al. (2011) and Nara et al. (2011) determined the permeability of the flow interior ranging from 10 -20 to 10 -19 m 2 . According to McGrail et al. (2006) and Zakharova et al. (2012), the presence of fractures or joints generally does not significantly contribute to the increase in basalt permeability. Figure 10 summarizes the statistics (min, max, mean, and standard deviation) of porosity and permeability estimations produced, respectively, by empirical Equations 1 and 2 for the different basalt hydrofacies.

CONCLUSIONS
The wireline logs presented in this study strongly suggest that basalts are best described as multilayered aquifers, in which permeable zones correspond to flow tops (vesicular, breccia, and scoria), that alternate with low permeability intervals represented by flow interiors (massive, columnar, and entablature basalt). Additionally, we recognized the presence of permeable intervals associated with weathered zones and clay lenses produced by intensive basalt weathering. Our analysis reinforces the use of well logs to identify the high and low permeability intervals of basalt flows. The analysis of compiled experimental data allowed us to produce an empirical equation capable to estimate porosity based on acoustic velocity logs. Our analysis also indicates that permeability may be successfully estimated from porosity using the Kozeny-Carman model. Thus, wireline logs provide a reasonable estimation of the hydraulic properties of basaltic aquifers.