Groundwater flow, heat flux, and permeability inferred from borehole temperature profile in Izu-Oshima volcano, Japan – Implications for subsurface hydrothermal system

Groundwater flow velocity as well as conductive and convective heat fluxes were estimated using temperature profile data from a 1000 m–deep borehole in the central part of the Izu-Oshima volcano, Japan. Two depth intervals below groundwater level with upward groundwater flow patterns were examined assuming a one-dimensional vertical steady flow. The groundwater velocity and total heat flux were estimated to be 5.0–5.4×10 -10 m/s and 0.54–0.59 W/m 2 , respectively, for the basement layer (Formation 4). For the shallower layer (Formation 2), both the upward velocity and heat flux were higher, indicating greater contributions of convective mass and heat transfer compared to those in the deeper layer. Furthermore, assuming that the upward flow was buoyancy-driven, vertical permeabilities of 2.8–5.1×10 -15 and 1.7–3.1×10 -16 m 2 , respectively, were estimated for Formations 2 and 4. The temperature patterns of the lava-dominant region (Formation 3), sandwiched between Formations 2 and 4, suggested the occurrence of lateral cooler groundwater inflow in fractures. These results were used for understanding a hydrothermal system beneath the volcano. The total heat flux estimated for Formation 4 (0.54–0.59 W/m 2 ) was nearly three times higher than the conductive heat flux in the northwestern coastal area, suggesting a higher heat supply below the central part of the volcano. A hydrothermal free convection system was inferred in Formations 2 and 3. In Formation 2, buoyancy-driven upward flow was enhanced because of the heat below and the higher permeability. Cooler groundwater was laterally supplied in lava fractures in Formation 3 to compensate for the mass loss by the upward flow at the bottom of Formation 2.


Introduction
Various volcanic phenomena are caused by interactions of magma and/or volatiles and their thermal energy with the background host medium and groundwater. Thus, it is important to understand the distributions and flow regimes of subsurface fluids and the physical properties of the host medium. In situ borehole temperature profile data have been widely used to derive information on subsurface heat and groundwater fluxes.
Particularly from non-linear temperature profiles, groundwater velocity and related convective heat flux as well as conductive one can be extracted (e.g. Anderson 2005).
Furthermore, if the driving force of the groundwater flow is known, the host medium permeability, which is crucial to constraining groundwater flow and heat transfer regimes, can be estimated. Such information is important for understanding hydrothermal systems below volcanoes (e.g. Lowell 1991).
Izu-Oshima is a basaltic island volcano with an elevation of 764 m in Japan ( Figure 1a and 1b). A low-elevation water table occurs on the island: around sea level near the coast and 36 meters above sea level (m.a.s.l.) in the central part of the volcano (at the borehole used in this study; Figure 1a). Approximately above sea level, forced 5 convection of recharged meteoric water through infiltration into the unsaturated layer and discharge into the sea via a hydraulic gradient after reaching the water table are expected. Furthermore, near the coast, saline water lies beneath freshwater, suggesting invasion of seawater into the volcanic edifice. Accordingly, the permeability of the host medium has also been estimated (Onizawa et al. 2009 and references therein). Contrary, to extract information associated with volcanic activity, deeper regions below sea level need to be studied, but little information has been obtained from in situ observations. Honda et al. (1979) estimated the conductive terrestrial heat flux using the borehole temperature profile of the basement layer (500-700 meters below sea level, m.b.s.l) in the northwestern coastal area (Figure 1a). These authors found a heat flux of 0.20 W/m 2 from a temperature gradient of 0.221 K/m and average thermal conductivity of 0.917 W/(m K) in rock samples from the borehole.
The 1000 m-deep borehole mentioned above was drilled in the western part of the caldera of the Izu-Oshima volcano by the Earthquake Research Institute, University of Tokyo. This borehole was located approximately 1 km west of the central cone, Mt.
Mihara, and penetrated the central part of the volcano (Figure 1a). The temperature 6 profile below the groundwater level (GWL) showed a high gradient and upward convex patterns, implying upward groundwater flow. Such information could aid the extraction of groundwater and heat fluxes and understanding of the hydrothermal system associated with volcanic activity.
In this study, vertical groundwater velocities and conducive and convective heat fluxes in the central part of the Izu-Oshima volcano were investigated using temperature profile data below the GWL. Furthermore, the permeability of the strata was determined by assuming that upward flow is buoyancy-driven. Based on the results, a hydrothermal system below the volcano was proposed.

Borehole Data
The wellhead altitude of the 1000 m-deep borehole is 547 m, and the bottom depth is 453 m.b.s.l.

Geological Formations
Geological formations at the borehole have been described by Nakada et al. (1999). In this study, the following formation numbers were used (Figure 1c-1e): 7 Formation 1 (Ground Surface to 4.5 m.a.s.l.) Approximately above sea level, the volcanic edifice mainly consists of the alternation of unconsolidated scoriae, volcanic ashes and lavas, which are subaerial deposits of the Oshima volcano.
Formation 2 (4.5 m.a.s.l. to 142 m.b.s.l.) Tuff breccia dominates the volcanic edifice below sea level in Formation 2, but the edifice is still composed of erupted materials from the Oshima volcano. This lithofacies reflects the growth of the volcanic edifice below and/or around sea level. Lavas dominate Formation 3. Nakada et al. (1999) unified this lava-dominant region with the tuff breccia-dominant region (Formation 2) as a geological unit. However, the bulk densities and temperature profiles of Formation 3 deviate from those of Formation 2. Since such differences were important to this study, we treated these formations as distinct.
8 Formation 4  Formation 4 consists of subaqueous sedimentary rocks, including clay minerals, such as montmorillonite and sericite. This formation formed a basement underlying the volcano and is correlated with the basement in the northwestern coastal area, where the terrestrial heat flux was estimated by Honda et al. (1979).

Density Profile
Variations in bulk density with changing depth were estimated via a borehole gravimeter measurement (BHGM) at intervals of approximately 10 m (Figure 1c). The BHGM density profile correlates with the formations and lithofacies described above.
In Formations 1-3, depth ranges coinciding with massive rocks, such as lavas and dikes, had higher densities (maximum of 2530 kg/m 3 ), whereas pyroclastic rocks had lower densities (minimum of 1415 kg/m 3 ). Such density variations were mainly attributed to variations in porosity. In Formation 4, the density variations (with a minimum and maximum of 2053 and 2220 kg/m 3 , respectively) were small compared to those in the upper formations. The temperature gradient seemed to be affected by lithology: At depths where massive lavas occurred, the temperature gradient became smaller compared to that at depths where clastic rocks occurred ( Figure 1e). The thermal conductivity of rocks was 10 generally greater than that of pore fluids, and the bulk thermal conductivity subsequently increased with decreasing porosity. Therefore, the effect of porosity should be corrected to quantify the conductive heat flux and, if present, the convective heat flux caused by groundwater flow.
Hereafter, we focused on the temperature profile below the GWL (Formations 2 to 4) to extract heat and groundwater fluxes originating from volcanic activity.

Water Properties
Two groundwater types were expected below the Izu-Oshima volcano: freshwater of meteoric origin and saline water from seawater invasion. However, the depth distributions of these waters around the borehole were uncertain. Thus, properties of these groundwater types were individually examined. The water properties were calculated with equations of state for conditions determined using the observed temperatures and pressures calculated assuming a hydrostatic state (Pritchett 1995).

Bulk Thermal Conductivity
As for the thermal conductivity, use of in situ information from the target borehole is an ideal manner. Porosity dependence is also required to correct the effect of the lithological differences. However, the only available data of Izu-Oshima were from the basement in the northwestern coastal area (Honda et al. 1979), and porosity dependence has not been investigated. Therefore, we predicted the bulk thermal conductivity for each depth interval via an empirical porosity-thermal conductivity relationship. Using the bulk density ̅ inferred from BHGM (Figure 1c), the porosity for each depth was estimated as follows: where and are the densities of the rock matrix and pore water, respectively.
Furthermore, was assumed to be 3000 kg/m 3 considering the basalts, and was calculated using the equations of state. The estimated porosity profile for freshwater is shown in Figure 2a, and the results were almost identical for saline water.
12 Robertson (1988) developed an empirical relationship for the bulk thermal conductivity of a mixture of rock and pore fluid ̅ to porosity and high-conducting mineral content (%): where and are intercept values of ̅ at 100 and 0% of the porosity, respectively, and is the correction coefficient for the mineral content. These parameters were determined for each rock type and pore fluid. For water-saturated mafic igneous rocks, was 0.753 W/(m K), and was 1.841 W/(m K). The mineral was olivine, with an S value of 0.025 W/(m K). Oshima basalts typically contain less than 3% olivine and pyroxene phenocrysts (Kawanabe 1991). Thus, in this study, 0 and 5% of the olivine content were examined to include actual conditions within these bounds. Furthermore, the constant value of 0.917 W/(m K) reported by Honda et al. (1979) was considered. Figure 2b shows the profiles of the three types of bulk thermal conductivity.

Conductive Heat Flux
The conductive heat flux ( ) (positive upward) is expressed as follows: where and are the temperature and depth (positive downward), respectively. By multiplying the calculated bulk thermal conductivity ( Figure 2b) and observed temperature gradient (Figure 2c), the depth-dependent conductive heat flux could be calculated ( Figure 2d). If the heat transfer within the target depth interval occurred only via conduction without convection, the conductive heat fluxes shown in Figure 2d should be constant regardless of depth. However, even when the effect of porosity was corrected by Robertson's model, the depth variations in the conductive heat fluxes remained ( Figure 2c). This indicated that heat transfer was not only conductive and suggested the presence of convective heat transfer by groundwater flow. Therefore, we considered the temperature profile data in light of the groundwater flow.

Vertical Groundwater Velocity and Conductive and Convective Heat Fluxes
The temperature profiles for Formations 2 and 4 were considered in light of the vertical flow, because both showed upward convex patterns (Figure 1d). The temperature and its gradient in Formation 3 largely deviated, making it difficult to apply the simple vertical 14 flow model. A qualitative consideration for this depth interval will be made in later section. Bredehoeft and Papadopulos (1965) first developed a method for extracting the velocity of steady vertical groundwater flow from temperature profile data. For a uniform thermal conductive medium, the governing equation is as follows: where ( ) and are the heat capacity and groundwater velocity, respectively.
While the depth is downward positive, the velocity and heat fluxes were set to be upward positive, since all the parameters considered in this study were upward. The solution of eq. (4) in a non-dimensional form is as follows: where , , and are the top and bottom temperatures and total thickness of the target depth interval, respectively, and is a non-dimensional parameter coinciding with the Péclet number, which represents the ratio of convective to conductive heat flux within the target depth interval (Clauser and Villinger 1990).
Since the depth variation of the thermal conductivity was obtained, the variation should be reflected in the analysis. Shan and Bodvarsson (2004) developed an equation for a one-dimensional vertical steady flow in a multilayered thermal conductive medium. The solution was identical to that for the uniform thermal conductivity model (eqs. (5) and (6)) by replacing the uniform bulk thermal conductivity ̅ in eq. (6) with an effective value ̅ ( ). This could be expressed by a harmonic mean of the depthdependent bulk thermal conductivity ̅ weighted by the layer thickness : where is the layer number, and is the number of layers from the top to the target depth. Finally, the vertical groundwater velocity was estimated using the temperature data.
The convective heat flux ( ) was described using the estimated velocity: depending on the choice of an arbitrary reference temperature . In this study, in Formations 2 and 4 was defined as 15 °C assuming the atmospheric temperature to detect the heat supply associated with volcanic activity. The total heat flux , 16 represented by the sum of the conductive and convective heat fluxes (eqs. (3) and (8)), should be constant within a target depth interval because of the required energy conservation: We examined six combinations of the two water types and three bulk thermal conductivity models mentioned above. Regarding depth intervals of the analyses, the entire interval was used for Formation 4. In contrast, it was difficult to specify the interval for Formation 2, because the temperature profile may have been affected by recharged water in the shallower part and by localized deviations of Formation 3 in the deeper part. In this analysis, the upper boundary was set at 28 m.b.s.l., where the temperature gradient peaked, and three depths were tested for the lower boundary.
The calculation conditions and obtained best-fit parameters are summarized in the additional file (Tables S1-S4). All the results indicated upward groundwater and heat fluxes. Deviations in the estimated parameters under the six conditions were primarily due to differences in the thermal conductivity models not in the water properties. Since the results were generally similar in all cases, only the cases for 17 freshwater and 0% olivine content in Robertson's model and the mid-depth interval of Formation 2 are depicted in Figure 3 as an example. The calculated temperatures fit the observed ones well (Figure 3a), and the root mean square of the residuals was less than 0.4 °C in all calculations. The groundwater velocity and total heat flux of Formation 4 were estimated to be 5.0-5.4×10 -10 m/s (16-17 mm/year) and 0.54-0.59 W/m 2 (at a reference temperature of 15 °C), respectively (Figure 3b and 3c). Although the temperature in Formation 2 was lower than that in Formation 4, the convective and total heat fluxes in the former were higher than those in the latter because of the higher groundwater velocity. This indicated that the convective mass and heat transfer were more significant in Formation 2.

Vertical Permeability
Upward groundwater velocity was investigated using the temperature profile. If the driving force gradient of the upward flow was known, the vertical permeability of the host medium could be estimated via Darcy's law. The most plausible cause of upward flow was the buoyant force caused by density variations in the groundwater, since the 18 high temperature gradient was observed. By assuming that the upward flow was buoyancy-driven, Darcy's law could be written as follows: where is the vertical permeability, is the gravitational acceleration, and , , and are the viscosity and densities at the top and bottom groundwater depths, respectively. Because the viscosity variations within the target intervals were large owing to the temperature variations, permeability was calculated for the minimum, maximum, and average viscosities (Tables S1-S4). The vertical permeability of Formation 4 was estimated to be 1.7-3.1×10 -16 m 2 . For Formation 2, the estimated vertical permeability was 2.8-5.1×10 -15 m 2 for all three depth interval analyses; these values were 9-30 times higher than those for Formation 4 ( Figure 4).

Permeability
Permeability in the Izu-Oshima volcano has been inferred by previous research (Figure 4). Pumping tests revealed a permeability of 6.9×10 -11 -6.3×10 -9 m 2 for the southeastern foot area (Ministry of Agriculture, Forestry and Fisheries of Japan 1980Japan , 1986, and a permeability of 5.7-7.8×10 -9 m 2 was estimated from the tidal responses of the GWL in the northwestern area (Koizumi et al. 1998).

20
The vertical permeability of formations below the GWL (Formations 2 and 4) were investigated in this study (Figure 4). The horizontal permeability was higher than the vertical permeability. This anisotropy was primarily attributed to horizontally stratified deposits. Furthermore, the vertical permeability decreased with increasing depth, mainly because of the differences in stratum lithology. Unconsolidated alternating layers of coarse scoriae and volcanic ashes in Formation 1 resulted in the highest permeability, whereas decreasing permeability was expected in the tuff breccias of Formation 2, even though both formations are composed of erupted materials of the Oshima volcano. Formation 4 had the lowest permeability because this basement sedimentary rock included fine materials, such as clay minerals and silts.

Qualitative Consideration of Groundwater Flow in Formation 3
In the lava-dominant Formation 3, the temperature and its gradient deviated from the trends in the upper and lower formations (Figure 1d and 1e). Three localized deviations were identified based on the temperature gradient. In particular, at 223 m.b.s.l., the temperature reached a local minimum (causing a negative temperature gradient) in the background increasing trend. Such temperature minima (and maxima) could not be explained by one-dimensional vertical groundwater flow and required a flow with a horizontal component (Lu and Ge 1996).
Particularly, "localized" temperature deviations were considered in light of groundwater flow in fractures. For a temperature minimum in the profile to be reached, inflow groundwater temperature should be lower than the background temperature (Ge 1998). An unpublished report of borehole drilling identified numerous fractures in the lavas of Formation 3. Considering the observational and theoretical facts above, lateral cooler groundwater inflow occurred in lava fractures in Formation 3.

Implications for Hydrothermal System below Volcano
The analytical results and inferences of groundwater flow, heat flux, and permeability obtained under the simple assumptions in this study were similar to a hydrothermal system with free convection of groundwater in porous media (e.g. Elder 1965).

Conclusions
The groundwater flow, heat flux, and permeability below the GWL were inferred using temperature profile data from the central part of the Izu-Oshima volcano. In light of previous studies, this study found that the vertical permeability was lower than the horizontal permeability and was also lower for deeper formations. In the basement, although the heat transfer was mostly conductive due to the low permeability, more heat was supplied than in the coastal areas. Geothermal free convection developed in relatively permeable formations above the basement. Buoyancy-driven upward flow occurred in the tuff breccia formation because of heating below. In the lava-dominant formation sandwiched between the tuff breccia and basement, cooler groundwater was laterally supplied within lava fractures.

Ethics approval and consent to participate
Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
The data used in this study belong to the Earthquake Research Institute, University of Tokyo. Contact the institute for any request regarding the data (http://www.eri.utokyo.ac.jp/en/).

Additional file
Table S1 Calculation conditions and estimated best fit parameters for Formation 2 (long-depth interval).