Development of Site-Scale Conceptual Model Using Integrated Borehole Methods: Systematic Approach for Hydraulic and Geometric Evaluation

: Understanding the physical ﬂow mechanisms in aquifer systems is essential in effectively protecting groundwater resources and preserving subsurface environments from a wide range of contaminants. A conceptual model is a simpliﬁed representation of a groundwater system and gaining knowledge about the geological features and parameters controlling the ﬂow and transport processes is a crucial ﬁrst step towards properly constructing a site-scale conceptual model. In this study, we present a multi-step workﬂow that involves integrated borehole techniques to gain information concerning groundwater ﬂow. Measurements from core-scale to ﬁeld-scale enable us to better build a subsurface geological structure divided into the unconsolidated layer and the fractured bedrock. In addition, neutron logging and mercury injection capillary pressure techniques allow for the development of vertical porosity distribution in the alluvial layer. For fracture characterization, the fracture geometry is delineated using a series of borehole imaging techniques and single-hole tests to differentiate the individual permeable fractures from other hydraulically inactive fractures. Combining the hydraulic and geometric evaluations, the presence of large-scale connective fracture networks is identiﬁed. Our high-resolution three-dimensional (3D) site-scale conceptual model is expected to contribute to improving the reliability and availability of numerical groundwater models.


Introduction
Groundwater shortages and the contamination of various types on various scales pose a significant threat to human health, ecosystems, and the environment. A conceptual model is a simplified representation of a groundwater system that is based on geological, geophysical, hydrological, and hydrogeochemical information [1]. Hydrogeological conceptual models have played an important role in managing groundwater resources, assessing contamination sites, and establishing remediation plans. Moreover, conceptualization is considered an important factor that can affect groundwater modeling results [2]. Therefore, there is a need for a systematic approach to building high-resolution 3D conceptual models that can help enhance the accuracy of numerical groundwater modeling.
Gaining knowledge associated with the geological features and parameters controlling flow and transport processes is a crucial first step toward properly constructing a hydrogeological conceptual model. Groundwater flow and contaminant transport depend on the physical properties of the rock and fluid or the geological characteristics of whether it is alluvial or consolidated bedrock. Although the mechanisms of groundwater flow through porous media are well-documented, many uncertainties still exist regarding groundwater flow and solute transport in fractured rocks [3,4]. In rock formations, fractures or secondary porosity are usually developed, considering that newly-formed rocks generally have very low porosity and permeability [5]. Nevertheless, fractures may act as conduits or barriers to groundwater flow, depending on their geometric and physical properties [6]. The key issue is, therefore, to provide more detailed insights into fracture networks in order to characterize their role in groundwater flow.
Site characterization often requires the assessment of a large number of fractures and the identification of preferential fluid flow-paths. Many in situ techniques including core inspection, various geophysical well loggings, and hydraulic tests have been used to identify permeable fractures [7][8][9][10]. Using these techniques, the hydrogeological properties of fractured formations are often determined by the discrete features of fractures and faults [11]. More specifically, individual fracture properties (e.g., aperture, length, and orientation), as well as properties of the fracture network (e.g., fracture density and connectivity), influence the hydraulic conductivity of crystalline fractured aquifers [12]. However, difficulties still exist when attempting to draw the geometry of fractures and their interconnections due to the random characteristics of their properties. For this reason, the comprehensive mapping of fracture networks typically takes significant time and effort [13].
High-resolution geophysical data can provide information at spatial scales and locations, and many previous studies have attempted to confirm the relationship between geophysical and hydrological parameters [14][15][16][17]. The temperature and electrical conductivity of groundwater are good indicators of groundwater flow as well as water quality [18,19]. In particular, it is possible to increase the sensitivity of identifying permeable intervals by using a differential-temperature log, which is a technique for recording the rate of temperature change with depth [20]. Borehole image logs record a continuous, magnetically oriented image of the borehole wall, making it possible to directly observe lithological changes and existing fractures in the surrounding rock [21,22]. Additionally, these down-hole imaging systems can construct effective fracture networks by measuring fracture orientation, dip, and aperture. However, as fracture density or aperture does not always dominate hydrogeological flow [23], fracture geometric data should be coupled with hydraulic data measured by a flowmeter or via a tracer test for the qualitative conceptualization of flow and transport [24]. Packer tests using inflatable packers are widely used in fractured rock hydrology to describe borehole conditions [25,26]. Doughty et al. [27] applied multi-rate flowing fluid electric conductivity logging to estimate the flow rate, transmissivity, and hydraulic head of individual fractures or high permeability zones. Dimmen et al. [28], referring to a wide range of the published literature, emphasize that fluid flow is primarily controlled by structure type, geometry, connectivity, kinematics, and interactions within the fracture network.
Despite numerous efforts to guide site characterization, the description of groundwater flow is still a challenging problem because it is 'site specific' and 'structurally heterogeneous'. Moreover, the type and extent of subsurface information required vary depending on the size and hydrogeological condition of the field site. Therefore, the objective of this study is to characterize the hydraulic properties and groundwater flow processes and to construct a high-resolution 3D conceptual model of groundwater flow and transport based on integrated borehole methods. For effective large-scale site description, we focus on characterizing the transmissive fractures and conduit connectivity and then highlight preferential flow-paths from the hydraulic and geometric characterization of the testbed site.

Geographical Setting
As part of the research project 'Subsurface Environment Management (SEM)', a dedicated Deokso testbed was established in Namyangju-si, South Korea to apply field testing techniques and analysis methods for characterizing subsurface flow and transport properties. The terrain of the region is characterized by low slopes, with mountains and wide rivers in the eastern and western regions. The center is dominated by a wide alluvial plain covered with thick weathered soils. To systematically apply borehole techniques, five boreholes were drilled, which were laid out in a square-shaped spacing pattern with approximately 40 to 80 m between boreholes ( Figure 1).

Geographical Setting
As part of the research project 'Subsurface Environment Management (SEM)', a dedicated Deokso testbed was established in Namyangju-si, South Korea to apply field testing techniques and analysis methods for characterizing subsurface flow and transport properties. The terrain of the region is characterized by low slopes, with mountains and wide rivers in the eastern and western regions. The center is dominated by a wide alluvial plain covered with thick weathered soils. To systematically apply borehole techniques, five boreholes were drilled, which were laid out in a square-shaped spacing pattern with approximately 40 to 80 m between boreholes ( Figure 1).
All boreholes were drilled at diameters of three or six inches and were up to approximately 100 m in depth. They were cased with steel to prevent collapse and then grouted with cement to about 50 m depth from the ground surface. The details of the boreholes are given in Table 1. All boreholes were fully cored. By interpreting these cores, the subsurface lithology was identified to be mostly fill, colluvium, residual soil, granite, and gneiss.   All boreholes were drilled at diameters of three or six inches and were up to approximately 100 m in depth. They were cased with steel to prevent collapse and then grouted with cement to about 50 m depth from the ground surface. The details of the boreholes are given in Table 1. All boreholes were fully cored. By interpreting these cores, the subsurface lithology was identified to be mostly fill, colluvium, residual soil, granite, and gneiss.

Methods
We present a multi-step workflow that involves various borehole techniques to gain information on the hydrogeological conditions of the testbed site ( Figure 2). A suite of integrated geophysical logging techniques, including gamma, spectral gamma, resistivity, velocity, fluid temperature, and electrical conductivity, enabled us to collect in situ data over the entire depth of each borehole ( Table 2). In addition to the field-scale data, seven cores representing each stratum of the testbed were sampled and analyzed to obtain core-scale physical properties such as permeability, porosity, density, capillary pressure, and velocity (Table 3)

Methods
We present a multi-step workflow that involves various borehole techniques to gain information on the hydrogeological conditions of the testbed site ( Figure 2). A suite of integrated geophysical logging techniques, including gamma, spectral gamma, resistivity, velocity, fluid temperature, and electrical conductivity, enabled us to collect in situ data over the entire depth of each borehole ( Table 2). In addition to the field-scale data, seven cores representing each stratum of the testbed were sampled and analyzed to obtain corescale physical properties such as permeability, porosity, density, capillary pressure, and velocity (Table 3).
Various borehole techniques were used to collect the structural information of the testbed, which was then constructed by geotechnical and geological approaches. For the geotechnical structure, we mainly used methods that can determine the hardness of the stratum, such as the sonic and density logs, while the geological strata were classified using the natural and spectral gamma ray logs that reflect the constituents of each stratum. In addition, obtained vertical log profiles were combined with core log interpretation which results in better definitions of both geological structure and lithological boundaries.   Various borehole techniques were used to collect the structural information of the testbed, which was then constructed by geotechnical and geological approaches. For the geotechnical structure, we mainly used methods that can determine the hardness of the stratum, such as the sonic and density logs, while the geological strata were classified using the natural and spectral gamma ray logs that reflect the constituents of each stratum. In addition, obtained vertical log profiles were combined with core log interpretation which results in better definitions of both geological structure and lithological boundaries.  To facilitate strategic planning, the approach to site characterization is divided into two categories: unconsolidated layer and fractured bedrock. Neutron log profiles were used to determine the in situ porosity governing groundwater flow in the unconsolidated layers, which were complemented by core porosity measured by a mercury injection porosimeter. The thermal neutron log, which measures the hydrogen index (HI) response, is expressed in apparent water-filled porosity units by assuming a constant matrix lithology [29][30][31] and is a very effective method for measuring in situ porosity in cased holes. Generally, the dual neutron log is widely used for porosity evaluation, and the ratio of neutron count rates from near and far detectors allows for the calculation of the porosity of the formation. This dual neutron log has the advantage of minimizing the influence of various borehole characteristics including borehole size, tool eccentricity, and borehole slant [32,33]. However, as steel casings with an inside diameter of 50 mm were installed for the purpose of hole stability in the boreholes, we used a single neutron sonde with a diameter of 41.3 mm. A single neutron log does not provide a calibration curve that converts the measured count rate into the formation porosity. Therefore, a calibration test was conducted to compensate for various problems including porosity conversion, the half-life of the neutron source, and the efficiency degradation of the detector. Single neutron log data were obtained from large-size blocks, including water, granite, and limestone in the radioactivity calibration facility. The thermal neutron count rates were measured in the borehole where the porosity is known from a dual neutron sonde [34][35][36]. These procedures helped establish the relationship between the count rate and the neutron porosity ( Figure 3). 022, 14, x FOR PEER REVIEW 6 of 20 Figure 3. The relationship between the count rate and the neutron porosity was derived from calibration tests at the radioactivity calibration facility (water, granite, and limestone blocks) and two test boreholes.
Understanding the distribution and characteristics of fractures intersecting boreholes is the beginning of fluid flow identification in complex fractured bedrock. A series of borehole imaging techniques, including an optical televiewer, caliper logging, and X-ray CT, allowed us to delineate the fracture geometry of the testbed with high resolution ( Figure  2). Fractures along the borehole were classified into six types: major open fracture, minor open fracture, partially open fracture, filled fracture, bedding/foliation, and broken zone. Additionally, each structural planar feature was statistically analyzed with respect to its azimuth, dip angle, and aperture. Vertical and spatial variations of groundwater level, temperature, and electrical conductivity were delineated from temperature and electrical conductivity logs. Temperature and electrical conductivity were measured at a speed of 3 m/min in each borehole at the ambient condition in equilibrium with the surrounding fluids. An attempt was made to distinguish fractures that contributed to groundwater flow under ambient conditions. Based on the information on fracture geometry and natural groundwater flow, a single-hole test was conducted to determine permeable fractures in each borehole ( Figure  4a). Vertical monitoring of temperature and electrical conductivity was conducted during a short-term pumping test. The pump used for the single-hole hydrogeophysical test is the variable frequency drive submersible pump, Rediflo2 (Grundfos), which is known to allow constant control of the flow rate by adjusting the frequency [37]. To effectively check the monitoring results within a short time, two pumps were placed at the lower part of the casing and the temperature and electrical conductivity of the fluid were monitored over the entire depth during constant pumping at a flow rate of 14~15 L/min (340 Hz).
After identifying hydraulically active fractures in each borehole, the multi-hole hydrogeophysical test was set to determine the connectivity of the fractured bedrock at the testbed site. The distance between the boreholes was about 40 m from borehole BH5-1, Figure 3. The relationship between the count rate and the neutron porosity was derived from calibration tests at the radioactivity calibration facility (water, granite, and limestone blocks) and two test boreholes.
Understanding the distribution and characteristics of fractures intersecting boreholes is the beginning of fluid flow identification in complex fractured bedrock. A series of borehole imaging techniques, including an optical televiewer, caliper logging, and X-ray CT, allowed us to delineate the fracture geometry of the testbed with high resolution (Figure 2). Fractures along the borehole were classified into six types: major open fracture, minor open fracture, partially open fracture, filled fracture, bedding/foliation, and broken zone. Additionally, each structural planar feature was statistically analyzed with respect to its azimuth, dip angle, and aperture. Vertical and spatial variations of groundwater level, temperature, and electrical conductivity were delineated from temperature and electrical conductivity logs. Temperature and electrical conductivity were measured at a speed of 3 m/min in each borehole at the ambient condition in equilibrium with the surrounding fluids. An attempt was made to distinguish fractures that contributed to groundwater flow under ambient conditions. Based on the information on fracture geometry and natural groundwater flow, a singlehole test was conducted to determine permeable fractures in each borehole (Figure 4a). Vertical monitoring of temperature and electrical conductivity was conducted during a short-term pumping test. The pump used for the single-hole hydrogeophysical test is the variable frequency drive submersible pump, Rediflo2 (Grundfos), which is known to allow constant control of the flow rate by adjusting the frequency [37]. To effectively check the monitoring results within a short time, two pumps were placed at the lower part of the casing and the temperature and electrical conductivity of the fluid were monitored over the entire depth during constant pumping at a flow rate of 14~15 L/min (340 Hz). which was located in the center (Figure 4b). Due to the long distance, a large capacity pumping rate was required to induce time-transient flow in the other monitoring wells. The submersible pump (Wilo PSB-L1012GP) was installed at the depth of the permeable zone in borehole BH5-1, where pumping was carried out for a total of 40 h at a pumping rate of 40 L/min. In the monitoring wells, the background vertical profiles of temperature and electrical conductivity were acquired before pumping, and changes were recorded by monitoring during and after pumping. For the characterization of groundwater flow within the boreholes, the vertical flow velocities can be measured by flowmeter loggings [38][39][40][41]. Since the groundwater flow rate is relatively slow at the testbed, a heat-pulse flowmeter [42] with a heat source was used to increase the sensitivity. The inter-borehole connectivity implemented through the hydraulic geophysical technique and the geometry (direction, inclination, and aperture) of the permeable fracture in each borehole were comprehensively analyzed in order to identify hydraulically active fractures.

Lithological Characterization
The geological distribution at the testbed can be identified with a natural gamma ray log. In this study, we tried to understand the geological connectivity of the high natural gamma ray zones and the detailed stratum characteristics through the spectral gamma ray log by measuring potassium (K), uranium (U), and thorium (Th), respectively. The intensity of the natural gamma rays at the testbed site showed an average of 137 API. The alluvial layer showed lower natural gamma rays than the bedrock layer, except for high gamma rays in the buried layer within 5 m of the surface. In BH5-1, the overall pattern was similar to that of BH2-1 and BH3-1 where high natural gamma rays were measured in some dyke sections and fractured zones. Table 4 summarizes the natural gamma ray After identifying hydraulically active fractures in each borehole, the multi-hole hydrogeophysical test was set to determine the connectivity of the fractured bedrock at the testbed site. The distance between the boreholes was about 40 m from borehole BH5-1, which was located in the center (Figure 4b). Due to the long distance, a large capacity pumping rate was required to induce time-transient flow in the other monitoring wells. The submersible pump (Wilo PSB-L1012GP) was installed at the depth of the permeable zone in borehole BH5-1, where pumping was carried out for a total of 40 h at a pumping rate of 40 L/min. In the monitoring wells, the background vertical profiles of temperature and electrical conductivity were acquired before pumping, and changes were recorded by monitoring during and after pumping.
For the characterization of groundwater flow within the boreholes, the vertical flow velocities can be measured by flowmeter loggings [38][39][40][41]. Since the groundwater flow rate is relatively slow at the testbed, a heat-pulse flowmeter [42] with a heat source was used to increase the sensitivity. The inter-borehole connectivity implemented through the hydraulic geophysical technique and the geometry (direction, inclination, and aperture) of the permeable fracture in each borehole were comprehensively analyzed in order to identify hydraulically active fractures.

Lithological Characterization
The geological distribution at the testbed can be identified with a natural gamma ray log. In this study, we tried to understand the geological connectivity of the high natural gamma ray zones and the detailed stratum characteristics through the spectral gamma ray log by measuring potassium (K), uranium (U), and thorium (Th), respectively. The intensity of the natural gamma rays at the testbed site showed an average of 137 API. The alluvial layer showed lower natural gamma rays than the bedrock layer, except for high gamma rays in the buried layer within 5 m of the surface. In BH5-1, the overall pattern was similar to that of BH2-1 and BH3-1 where high natural gamma rays were measured in some dyke sections and fractured zones. Table 4 summarizes the natural gamma ray results at each borehole. The testbed site is mainly composed of gneiss, and the natural gamma rays are generally higher in the west (BH1-1 and BH4-1) than in the east (BH2-1 and BH3-1). This indicates that geological characteristics between the west and east are different (Figures 1 and 5). The K, U, and Th values were calculated using the energy spectrum from spectral gamma ray logging (Figure 6a). As a result of analyzing each element in zone A, where the section showed the high natural gamma ray in the vertical profiles, it was found that the geological connectivity between the center area (BH5-1) and the west area (BH1-1 and BH4-1) was insufficient (Figure 6b). This indicates that there is a difference in the composition of the strata, unlike the strata structure interpreted only by natural gamma ray responses (Table 5). To determine the P-wave velocity of the strata, sonic loggings were performed and data were acquired with a sampling rate of 4 µs using a 15 kHz frequency transmitter and three receivers. P-wave velocities were obtained only in the open section of the borehole. Overall, BH1-1 and BH4-1 showed a high-velocity distribution, and there was a difference in the trend between the east (BH1-1 and BH4-1) and west (BH2-1, BH3-1, and BH5-1) of the testbed (Figure 7a). The characteristic difference between these eastern and western geological features are similar to the observations from the natural gamma ray analysis, supporting geological differences between the eastern and western regions. Additionally, it was confirmed that a low-velocity layer existed in BH2-1, BH3-1, and BH5-1 at greater depths (Figure 7b). These velocity profiles ensure a vertical distribution of weathered rocks, soft rocks, and hard rocks. different (Figures 1 and 5). The K, U, and Th values were calculated using the energy spectrum from spectral gamma ray logging (Figure 6a). As a result of analyzing each element in zone A, where the section showed the high natural gamma ray in the vertical profiles, it was found that the geological connectivity between the center area (BH5-1) and the west area (BH1-1 and BH4-1) was insufficient (Figure 6b). This indicates that there is a difference in the composition of the strata, unlike the strata structure interpreted only by natural gamma ray responses (Table 5).     To determine the P-wave velocity of the strata, sonic loggings were performed and data were acquired with a sampling rate of 4 µs using a 15 kHz frequency transmitter and three receivers. P-wave velocities were obtained only in the open section of the borehole. Overall, BH1-1 and BH4-1 showed a high-velocity distribution, and there was a difference in the trend between the east (BH1-1 and BH4-1) and west (BH2-1, BH3-1, and BH5-1) of the testbed (Figure 7a). The characteristic difference between these eastern and western geological features are similar to the observations from the natural gamma ray analysis, supporting geological differences between the eastern and western regions. Additionally, it was confirmed that a low-velocity layer existed in BH2-1, BH3-1, and BH5-1 at greater depths (Figure 7b). These velocity profiles ensure a vertical distribution of weathered rocks, soft rocks, and hard rocks. Analyzing the natural/spectral gamma rays and velocity, we reclassified the core logs and developed a 3D subsurface structure of the testbed (Figure 8a). The alluvial layer is widely distributed in the direction of BH1-1 and BH4-1, and the bedrock is distributed with a slope from west to east (Figure 8b,c). Analyzing the natural/spectral gamma rays and velocity, we reclassified the core logs and developed a 3D subsurface structure of the testbed (Figure 8a). The alluvial layer is widely distributed in the direction of BH1-1 and BH4-1, and the bedrock is distributed with a slope from west to east (Figure 8b,c). Analyzing the natural/spectral gamma rays and velocity, we reclassified the core logs and developed a 3D subsurface structure of the testbed (Figure 8a). The alluvial layer is widely distributed in the direction of BH1-1 and BH4-1, and the bedrock is distributed with a slope from west to east (Figure 8b,c).

Unconsolidated Layer Characterization
The porosity of the testbed was estimated using the empirical formula derived from the single neutron responses measured from the calibration pits. The measurement interval was 10 cm and the total porosity of the strata was calculated after removing the casing effect. The range of porosity in the alluvial layer was 20%-42% (Figure 9a), whereas the rock layer below 50 m had a porosity smaller than 10%. The porosity increase corresponded to the presence of the weathering zone and the fractured zone, for the most part (Figure 9b). The porosities and pore size distributions of the core samples (Table 3) were analyzed using a mercury injection porosimeter (Micromeritics Autopore system) for the purposes of measuring the representative porosity at the unsaturated zone and confirming the reliability of the neutron logging predicted porosity. The core porosity is the effective porosity from measuring only the volume of pores connected to each other in the core sample. Here, the relationship between the pore diameter and the pressure is calculated using Washburn's equation [43]. As a result of the core porosity measurement, it was confirmed that the porosity decreased with increasing depth (Table 6). By comparing the porosity obtained from the core test with the vertical porosity profile derived from the neutron logs, it was confirmed that the trend of the results fits well even if the value of the core test is the effective porosity and the neutron logs are taken as the total porosity (Figure 9c). Table 6. Core porosity measuring from the mercury injection capillary pressure method; the sampled depths and lithologic information in BH3-1 are noted in Table 3. sample. Here, the relationship between the pore diameter and the pressure is calculated using Washburn's equation [43]. As a result of the core porosity measurement, it was confirmed that the porosity decreased with increasing depth (Table 6). By comparing the porosity obtained from the core test with the vertical porosity profile derived from the neutron logs, it was confirmed that the trend of the results fits well even if the value of the core test is the effective porosity and the neutron logs are taken as the total porosity (Figure 9c).

Fractured Bedrock Characterization
From the porosity estimation, the bedrock at the testbed showed a negligible porosity of less than 1%, suggesting that fractures serve as the primary passage for groundwater flow. In order to evaluate the distribution and characteristics of the fractures, an optical televiewer was used. Then, the lithology identification, fracture distribution, fracture development, classification, and aperture size were analyzed at a high resolution. BH1-1, BH4-1, and BH5-1, located in the west of the testbed, showed a large number of discontinuities; the azimuth of each fracture zone was formed in the southeast direction with a dominant incline of 30-40 • (Figure 10a). The aperture size of the fracture was 10 mm or less with the presence of multiple large fractured zones (over 200 mm) ( Table 7). The caliper log, which measures the diameter of the borehole using three arms, was only able to measure large fractures. However, by analyzing the optical televiewer images, it was identified that open fractures where fluid movement can occur are distributed in BH1-1, 2-1, and 5-1. In addition, some partially filled fracture zones are mainly concentrated in BH3-1 and 4-1. The fracture distribution in the borehole was visualized after analyzing the optical images and caliper logs, and fracture density was identified with fracture height analysis (Figure 10b).

Groundwater Flow Characterization
Through the characterization of the fractured bedrock at the testbed site, the permeable fracture zones at each borehole were identified by integrating the fracture zone size from the caliper logs, the velocity attenuation from the sonic logs, and the inflection point from the differential-temperature logs ( Figure 11). This is an approximate result to quickly estimate groundwater flow using the borehole techniques under ambient flow conditions.
A single-hole test was conducted to differentiate individual permeable fractures from other hydraulically inactive fractures. After acquiring the background vertical profiles under ambient flow conditions, the groundwater level and temperature/electrical conductivity were measured at 30-min intervals during a pumping test. They were also monitored even after the pumping was stopped to examine the recovery phenomenon. As a result of monitoring, it was revealed that the change in electrical conductivity was more sensitive than the change in temperature. Figure 12a shows the monitoring results of the single-hole test in BH5-1, and the primary groundwater flow was identified at two points (54 m and 57.5 m). Both these permeable fracture zones ( 1 and 2 ) are major open fractures with an inclination of about 30-50 • to the SW azimuth (Figure 12b).
At the permeable fracture zones of BH5-1, electrical conductivity distinctively fluctuated during pumping activity. Based on this trend, these fracture zones were selected as target fracture zones for multi-hole hydraulic geophysical testing which is able to determine the connected fractures through boreholes in the testbed. During the pumping activity at BH5-1, only BH2-1 and BH3-1 showed a drop and recovery of groundwater level; such changes imply that both boreholes are hydraulically connected to BH5-1. On the other hand, no changes in groundwater levels were observed in BH1-1 and BH4-1 ( Figure 13).
Electrical conductivity did not change at BH1-1 and BH4-1, but distinct changes were detected at some depths (BH2-1: 50 m and 85 m and BH3-1: 65 m) (Figure 14a,b). As can be seen, pumping activity at BH5-1 led to groundwater inflow and outflow through connected fracture zones at both BH2-1 and BH3-1 and induced changes in both groundwater level and electrical conductivity. To further confirm the connectivity of the fracture zones, a flowmeter analysis was performed. The measurement depth was targeted at the upper and lower parts of the known fracture zones where connectivity was inferred through temperature/electrical conductivity monitoring. BH2-1 showed no flow under ambient conditions, but a down-flow (negative value) occurred during pumping in BH5-1 (Figure 14a). Similarly, BH3-1 showed no significant change under ambient conditions, but an up-flow occurred with a slight increase in the flow rate (Figure 14b). Note that no significant changes in the flow rate were noticed in BH1-1 or BH4-1, even under pumping conditions. 2022, 14, x FOR PEER REVIEW 21 of 31 Figure 11. Vertical profiles from temperature, differential temperature, caliper, and sonic logs in boreholes BH1-1, BH2-1, A single-hole test was conducted to differentiate individual permeable fractures from other hydraulically acquiring the background vertical profiles under ambient flow conditions, the groundwater level and tempera tivity were measured at 30-min intervals during a pumping test. They were also monitored even after the pu Figure 11. Vertical profiles from temperature, differential temperature, caliper, and sonic logs in boreholes BH1-1, BH2-1, BH3-1, and BH4-1.
x FOR PEER REVIEW 22 of 31 examine the recovery phenomenon. As a result of monitoring, it was revealed that the change in electrical conduct sensitive than the change in temperature. Figure 12a shows the monitoring results of the single-hole test in BH5-1, a groundwater flow was identified at two points (54 m and 57.5 m). Both these permeable fracture zones (① and ②) fractures with an inclination of about 30-50° to the SW azimuth (Figure 12b).
(a) (b) Electrical conductivity did not change at BH1-1 and BH4-1, but distinct changes were detected at some depths (BH2-1: 50 m and 85 m and BH3-1: 65 m) (Figure 14a,b). As can be seen, pumping activity at BH5-1 led to groundwater inflow and outflow through connected fracture zones at both BH2-1 and BH3-1 and induced changes in both groundwater level and electrical conductivity. To further confirm the connectivity of the fracture zones, a flowmeter analysis was performed. The measurement depth was targeted at the upper and lower parts of the known fracture zones where connectivity was inferred through temperature/electrical conductivity monitoring. BH2-1 showed no flow under ambient conditions, but a down-flow (negative value) occurred during pumping in BH5-1 ( Figure  14a). Similarly, BH3-1 showed no significant change under ambient conditions, but an upflow occurred with a slight increase in the flow rate (Figure 14b). Note that no significant changes in the flow rate were noticed in BH1-1 or BH4-1, even under pumping conditions. These connectivity results were combined with the geometrical information, such as  These connectivity results were combined with the geometrical information, such as the azimuth, dip direction, and aperture of individual fractures. It was thus confirmed that fractures at a depth of approximately 53 m in BH5-1 are geometrically connected to permeable fractures in BH2-1 and BH3-1 (Table 8). Consequently, a field-scale 3D hydraulic characterization model representing the testbed was constructed ( Figure 15). The model includes a porous unconsolidated layer and a fractured bedrock. The porosity in the unconsolidated layer was characterized by a neutron porosity that can be applied even in the cased boreholes. In addition, the procedure for establishing the fracture model allows us to identify the characteristics and distribution of the fractures and estimate the permeable fracture zone at each borehole as well as aid us in evaluating the geometrically and hydraulically connective fractures.    Through the construction of the 3D high-resolution site-scale conceptual model, the major conclusions are listed below. First of all, based on the analyses of the gamma-ray responses and the velocity distribution, the testbed shows geological heterogeneity, especially between the east and west sides (Figures 5-7). The hydraulic conductivity, groundwater level drop, and short-term pumping test highlight a distinct difference between the west area (BH1-1 and BH4-1) and the east area (BH2-1 and BH3-1) ( Table 9). Changes in the groundwater level during the pumping test at BH5-1 further support the existence of a hydraulic connection ( Figure 13). Additionally, from the borehole-to-borehole electrical resistivity tomography analyses conducted prior to the pumping test, it was possible to understand the overall development of the subsurface fracture zones. The high resistivity distribution between the center area (BH5-1) and the west area (BH1-1 and BH4-1) acts as a hydraulic barrier, and the low resistivity anomalies connected between the center area (depth of 50 m in BH5-1) and the east area (BH3-1 (Figure 16a) and BH2-1 (Figure 16b)) suggest the possible existence of connected fractures. Through the construction of the 3D high-resolution site-scale conceptual model, the major conclusions are listed below. First of all, based on the analyses of the gamma-ray responses and the velocity distribution, the testbed shows geological heterogeneity, especially between the east and west sides (Figures 5-7). The hydraulic conductivity, groundwater level drop, and short-term pumping test highlight a distinct difference between the west area (BH1-1 and BH4-1) and the east area (BH2-1 and BH3-1) ( Table 9). Changes in the groundwater level during the pumping test at BH5-1 further support the existence of a hydraulic connection ( Figure 13). Additionally, from the borehole-to-borehole electrical resistivity tomography analyses conducted prior to the pumping test, it was possible to understand the overall development of the subsurface fracture zones. The high resistivity distribution between the center area (BH5-1) and the west area (BH1-1 and BH4-1) acts as a hydraulic barrier, and the low resistivity anomalies connected between the center area (depth of 50 m in BH5-1) and the east area (BH3-1 (Figure 16a) and BH2-1 (Figure 16b)) suggest the possible existence of connected fractures.

Discussion
The properties of intrinsic rock and fluid can influence fluid flow. Additionally, geological structures such as fractures, faults, and deformation bands are also considered important factors affecting fluid flow characteristics [44,45]. Lods et al. [46] conducted various kinds of pumping experiments at two 15 m-apart boreholes to estimate complex hydrogeological features. In addition, Paillet et al. [47] identified the fracture connections in Figure 16. A 2D hole-to-hole cross-section from electrical resistivity tomography: (a) BH1-1, BH5-1, and BH3-1; (b) BH4-1, BH5-1, and BH2-1. Dotted red lines denote inferred connective fractures derived from the 3D conceptual model.

Discussion
The properties of intrinsic rock and fluid can influence fluid flow. Additionally, geological structures such as fractures, faults, and deformation bands are also considered important factors affecting fluid flow characteristics [44,45]. Lods et al. [46] conducted various kinds of pumping experiments at two 15 m-apart boreholes to estimate complex hydrogeological features. In addition, Paillet et al. [47] identified the fracture connections in two boreholes spaced 21 m apart in Melechov granite in the Bohemian-Moravian Highland, Czech Republic. Despite their efforts, these studies suggest that the assessment of hydraulic connectivity is very difficult and time-consuming, especially for a large-scale aquifer system.
Compared to previous studies, the size of the testbed in this study is large (80 × 80 m) with a distance of 40 m between boreholes. Nevertheless, by implementing multiple hydrogeologic and geophysical borehole logging techniques, we effectively developed a 3D high-resolution site-scale conceptual model. The major in situ physical properties were measured from the core-scale to the site-scale, and the hydraulic properties were evaluated. The porosity of the strata is typically measured in the laboratory with core samples collected at the field site, and thus, it cannot take into account the conditions of the strata influenced by in situ pressure. To overcome such limitation, neutron logging, density logging, and sonic logging was applied to estimate the in situ total porosity [48,49]. Meanwhile, for the characterization of the fractured rock, the geometric configuration of the fractures at the testbed was quantified via analyzing high-resolution images, and individual permeable fracture zones were verified through large-scale hydraulic tests and monitoring. In addition, the connectivity of groundwater flow was also demonstrated based on hydraulic and geometric evaluations. Even though subsurfaces are complex with a heterogeneous nature and existing fractures, it was revealed that the application of multiple hydrogeologic and geophysical borehole techniques aided the understanding of fracture connectivity and associated groundwater flow throughout the study.

Conclusions
A variety of borehole techniques have been devoted to constructing a 3D hydraulic geological conceptual model for understanding the behavior of subsurface groundwater. Although the testbed had many obstacles due to the geological heterogeneity, large scale, and a small number of boreholes, an attempt was made to construct a high-resolution 3D conceptual model using a systematic approach. The vertical profiles of the in situ physical properties and the laboratory core test results enabled us to effectively build a geological model, which is composed of both the unconsolidated layer and the fractured bedrock. The model focused on identifying the distribution of porosity and fractures as the main factors responsible for groundwater flow. The porosity was estimated using a single neutron probe available for small-sized boreholes, and calibration blocks were used for quality control. The porosity of the unsaturated zone where the neutron logging cannot be applied was supplemented using mercury injection capillary pressure techniques. For the hydraulic and geometric characterization of the fractures, the position, inclination, and direction of the discontinuities were statistically established through borehole techniques. Then, the groundwater flow model was proposed through hydraulic tests and monitoring.
Hydrogeological conceptual site models describe groundwater systems, which are considered a major source of uncertainty in groundwater flow and transport modeling [50]. We derived a site-scale conceptual model through a deterministic approach using in situ data, which is expected to contribute to improving the reliability of numerical modeling by providing actual input data for fluid-contaminant behavior modeling. In particular, even in the case of non-coring boreholes or boreholes with poor core recovery rates, our systematic approach can provide key information for numerical modeling.