Enhanced Characterization of Fractured Zones in Bedrock Using Hydraulic Tomography through Joint Inversion of Hydraulic Head and Flux Data

: Hydraulic tomography (HT) is a promising technique for high-resolution imaging of sub-surface heterogeneity, which addresses the limitations of traditional methods, such as borehole drilling and geophysical surveys. This study focuses on the application of HT to detect and characterize fractured zones in bedrock and addresses the gap in the understanding of the role of distributed ﬂ ux data in the joint inversion of hydraulic head and ﬂ ux data. By conducting synthetic injection tests and using sequential successive linear estimators for inversion, the study explores the effectiveness of combining limited head data with distributed temperature sensing (A-DTS)-derived ﬂ ux data. The ﬁ ndings highlight the fact that integrating ﬂ ux data signi ﬁ cantly enhances the accuracy of identifying fracture permeability characteristics, even when head data is sparse. This approach not only improves the resolution of hydraulic conductivity ﬁ elds but also o ﬀ ers a cost-e ﬀ ective strategy for practical ﬁ eld applications. The results underscore the potential of HT to enhance our understanding of groundwater ﬂ ow and contaminant transport in fractured media, which has important implications for carbon capture


Introduction
Fractures serve as preferential flow pathways, influencing fluid migration and solute transport in rocks.Understanding fracture structure and hydraulic properties is essential for predicting fluid migration and contaminant transport, which is crucial for a variety of applications, including carbon capture and storage (CCS), enhanced geothermal systems (EGS), and high-level radioactive waste disposal [1].Various methods, such as borehole drilling, geophysical surveys, and hydraulic experiments, have long been used to gather valuable information for fracture characterization [2,3].However, these methods often yield ambiguously averaged estimates that lack the resolution required for detailed analysis.Hydraulic tomography (HT), which has emerged as an efficient and powerful technique for high-resolution imaging of subsurface heterogeneity, addresses the limitations of traditional methods by utilizing the joint inversion of data from multiple pumping tests and multiple sources of data.
Various inversion methods, including successive linear estimators, sequential successive linear estimators, quasi-linear geostatistical approaches, and ensemble Kalman filtering, improve the robustness of HT in identifying hydraulic conductivity variations and offer a comprehensive approach to aquifer characterization with realistic and computationally efficient results [4][5][6][7].
Early applications of HT were primarily in unconsolidated deposits, but recent studies have extended its use to fractured rocks, and it shows promising results in studies based on numerical simulations, sandbox experiments, and field investigations [6,[8][9][10].Studies by Illman et al. [10] and Zha et al. [11] demonstrate the effectiveness of HT in mapping fracture patterns and fault zones in fractured granite, showcasing its potential for large-scale applications.The use of synthetic experiments by Hao et al. [12] and numerical simulations by Sharmeen et al. [13] further validate the capability of HT to map fractures and predict drawdown behavior accurately.Additionally, the integration of geophysical data and tracer tests, as shown by Sreeparvathy et al. [14], enhances the characterization of fracture-rock aquifers.Overall, HT proves to be a promising method to map fracture and fault distributions in fractured aquifers, aiding in better understanding groundwater flow and contaminant transport [15].
However, HT faces several challenges and limitations when applied to fractured aquifers.Accurately mapping the fractures and their hydraulic properties is difficult due to the inherent heterogeneity and complex connectivity of fractures [9].The resolution of HT results depends heavily on the density of pumping and monitoring locations and the quality of data collected [11,16].Moreover, the spatial representative elementary volume (REV) of fracture connectivity must be considered; if this REV does not exist, a dense observation network is required to achieve satisfactory results [17].Considering that both borehole drilling and the double-packer hydraulic testing system in fractured rocks are extremely costly, integration of other data sources, such as groundwater flux, temperature, and geophysical surveys, could be very important [13].
Flux data play a crucial role in hydraulic tomography by enhancing the estimation of hydraulic conductivity fields when used jointly with head data [18][19][20].The inclusion of flux measurements provides non-redundant information about aquifer heterogeneity, improving the resolution of estimated properties, especially in short-term pumping tests.Studies have shown that combining head and flux data leads to more accurate estimations of hydraulic conductivity, thus reducing uncertainty and enhancing the overall performance of hydraulic tomography [20,21].Additionally, flux data contribute significantly to parameter inversion, particularly when geological data are available, which makes their joint collection with head data highly recommended for comprehensive and effective hydraulic tomography surveys.
The monitoring of groundwater fluxes variations is challenging because traditional approaches, such as piezometers, tracer tests, or dilution tests, remain limited for continuous flux monitoring and only provide the spatially averaged value of fluxes, requiring repeating point measurements at different locations [22].However, active distributed temperature sensing (A-DTS), which has developed rapidly in recent years, provides highresolution spatial and temporal data on groundwater flux using heat as a tracer along fiber optic (FO) cables in various hydrogeological contexts, including bedrock aquifers, alluvial aquifers, and streambed sediments [22][23][24][25].By heating the FO cable and monitoring the temperature response, A-DTS can estimate groundwater flow rates using temperature as a tracer and provide a robust and versatile approach to monitoring and quantifying groundwater fluxes, which makes possible the joint inversion of groundwater head and flux, which can be used in practical applications for fracture characterization.
Despite these improvements of the joint inversion methods of hydraulic head and flux data, accurately characterizing fractured zones in complex geological settings remains challenging.How to take the greatest possible advantage of the flux data and maximize the non-redundant information of the flux data are two of the urgent problems to be solved in the practical application of HT.This study addresses the gap in understanding the role of flux data in joint inversion by synthetic experiments using a sequential successive linear estimator.The influence of flux data was explored and a strategy for its application in characterizing fractured zones in bedrock was developed.This paper is organized as follows: Section 2 describes the theoretical framework and methodology; Section 3 presents the numerical simulations results; Section 4 discusses the synthetic experiments and findings; and Section 5 concludes with a summary and implications for future research.

Materials and Methods
First, a two-dimensional model of the fractured zone was established, and synthetic injection tests were conducted to obtain datasets of groundwater heads and flow flux.Second, inversion models were constructed using the hydraulic tomography method, employing these head and flux datasets to estimate hydraulic conductivity.Finally, independent validation tests were conducted to compare the performance of inversion models using head, flux, and combined datasets.

Model Setup
A cross-sectional model of a fractured zone was established as a reference model using the software COMSOL Multiphysics.Synthetic injection tests were conducted to obtain water head and flux datasets.The model is conceptualized as a two-dimensional vertical cross-section of a specific site, with a depth of 500 m and a focused width of 100 m (Figure 1).To avoid boundary condition effects, the model extent is expanded to 500 m by 500 m.A horizontally oriented, irregularly shaped fractured zone passes through the model.The boundaries of the reference model are set as constant water head boundaries on the sides, while the top and bottom boundaries are zero-flux boundaries.There is a background groundwater flow from the left side (x = 0 m) to the right side (x = 500 m), with a hydraulic gradient of 0.01.The hydraulic conductance in the fractured zone and matrix is 5 × 10 −4 m/s and 5 × 10 −9 m/s, respectively.

Synthetic Injection Tests Used for Inversion and Validation
To obtain the water head and flux data required for joint inversion, three synthetic injection tests were conducted in the reference model.For each injection test, an injection point was set (specific injection points are shown in Table 1), with an injection rate of 10 kg/s (corresponding to a volumetric flow rate of 0.01 m 3 /s), and the injection duration was 10 days.From each injection test, water head response data from 15 monitoring points across 5 borehole profiles (including the injection point) and distributed flux profiles were obtained.
To obtain the necessary water head and flux data for joint inversion, we conducted three synthetic injection tests within the reference model.In each test, a single injection point was established (specific locations are detailed in Table 1), with a water injection rate of 10 kg/s (equivalent to a volumetric flow rate of 0.01 m 3 /s) over a period of 10 days.

Test ID Injection Point Head Dataset Flux Dataset
Well5-2 H is the water head; Vx is the horizontal flow rate; Vy is the vertical flow rate; n is the number of data.
Table 1 summarizes the synthetic injection tests conducted within the fracture zone reference model, specifying the injection points and the corresponding water head and flux datasets obtained.All data were collected at the end of the injection tests, ensuring that water levels and flow rates at the monitoring points had reached a quasi-steady state.The hydraulic monitoring points provided water head (H) data, while the flux monitoring points provided horizontal and vertical groundwater flow velocities (Vx and Vy, respectively).These datasets served as input sources for the hydraulic tomography inversion to identify fracture permeability characteristics.
Independent hydraulic tests in both the reference model and the estimated permeability distribution fields obtained from the inversion used different injection points from those in the inversion simulations to obtain water level (H) and flow rate (Vx and Vy) distributions for each test.Considering that the actual monitoring simulation range of the reference model is only 100 m × 500 m (expanded to 500 m × 500 m to avoid boundary condition effects), only the area of x = 150 m to 350 m and y = −500 m to 0 m is selected as the validation range (200 m × 500 m).The validation range is divided into a 5 m × 5 m grid (consistent with the VSAFT2 inversion model grid size), and the water levels and flow rates at each grid node (totaling 10,201 validation points) are compared and analyzed using validation scatter plots.
Independent hydraulic tests were conducted within both the reference model and the model with the estimated permeability fields derived from the inversion.These tests used different injection points at Well2-2 (inside the fractured zone), Well3-1 (in the matrix above the fractured zone), and Well4-3 (in the matrix beneath the fractured zone), which are different from those in the inversion simulations to obtain water head and flow.
Considering that the monitoring range of the reference model is 100 m × 500 m, validation was limited to x = 150 m to 350 m and y = −500 m to 0 m, covering an area of 200 m × 500 m.
Groundwater head and flux in the validation area, 10,201 validation points in total, were compared and analyzed.The performance was assessed using three indicators: the coefficient of determination (R 2 ), the mean absolute error (MAE), and the root mean square error (RMSE).

Inverse Modeling Approach
Hydraulic tomography of the synthetic injection test data from the reference model was conducted using the simultaneous successive linear estimator (SimSLE) algorithm [5].The SimSLE can estimate the K and SS distributions as well as their uncertainties by evaluating all the available pumping/injection test data simultaneously.A brief description of the SimSLE is provided below.
In the SimSLE, the hydraulic parameters of the spatial random field are generalized as the following expressions: In the above expressions, ln is the natural logarithm, E[ln K] and E[ln SS] are the mean values, and f and S are the perturbations of the natural logarithm values of the parameters.The instantaneous dynamic water head response caused by the pumping test can be expressed as the sum of the mean value and the perturbation, as follows: The SimSLE method requires starting calculations based on the prior values of the parameter field.The difference between the simulated and the observed water head at sampling points is used for the inversion of hydraulic parameter values in the next iteration.The parameter estimation field is iteratively generated as follows: In the equation, r  is the parameter estimation value at the r iteration, and the sub- script c represents specific conditions.The weight factor ω varies according to the difference between the observed and simulated values at each time node in each observation well.The superscript T represents the transpose matrix, * d is the vector of observed wa- ter head, and ( ) r d is the water head value calculated by the forward simulation at the r iteration.The weight factor ω, as an important calculation parameter, can be calculated by the following formula: In the above formula, r dd ε is the conditional covariance matrix between the second- ary datasets (i.e., head, flux magnitude, or flux vector) at the r iteration, and r dy ε is the cross-covariance matrix of the secondary datasets and calculation parameters at the r iteration.θ is the stability coefficient, which needs to be added to the diagonal matrix ( ε can be calculated by first-order approximation nu- merical methods [8] and can be written as follows: In the equation, Jd is the sensitivity matrix of the water head to the estimated parameter values.In the first iteration (i.e., r = 0), the covariance matrix εyy is composed of prior geological information, such as K and Ss.When the number of iterations is r ≥ 1, the conditional covariance matrix can be calculated by the following formula:

Inverse Model Setup
This study employs VSAFT2 (see https://tian.hwr.arizona.edu/downloadsaccessed on 10 June 2023) as the numerical code for the inversion model.A vertical profile inversion model is established in VSAFT2, maintaining consistency with the reference model i terms of the range, number, and locations of the water head and flux monitoring points.Given that the input dataset, water heads, and flux at the end of the simulated injection tests represents steady-state flow data, the selected inversion method is a steady-state flow inversion.
To evaluate the impact of different data types (head and flux) on the inversion, three different input datasets were used for hydraulic tomography:

•
Model_V uses only the flux dataset, utilizing data from all 105 flux monitoring points; • Model_H uses only the water head dataset, incorporating data from 15 monitoring points; • Model_H-V combines the water head and flux data using the results of Model_V as the prior parameter field and then integrating the water head data for joint inversion.

K Estimates from Flux Data and Model Validation
Figure 2 shows the inversion results of Model_V, based solely on the flux data.By comparing Figure 2a with Figure 1, it can be seen that the permeability distribution estimated from the flux data can largely depict the position of the high-permeability fractured zone.It also accurately reflects the significant morphological change, narrowing the width of the fractured zone at x = 275 m.However, the estimation of matrix permeability in Figure 2a is poor, with the matrix permeability differing significantly from that in the reference model by about three orders of magnitude.The errors in Figure 2a correspond to the uncertainty distribution shown in Figure 2b.The uncertainty is low at the positions where flux information was obtained in the boreholes, while the uncertainty is high in the surrounding matrix areas due to the lack of information.As shown in Figure 2c-e for Model_V, based solely on the 15 water heads monitoring data, water head H fits better for the injection test conducted at the monitoring point located in the high-permeability fractured zone (Well2-2) compared to the injection tests conducted in the matrix (the injection points at Well3-1 or Well4-3), with R 2 = 0.42892, MAE = 12.3 m, and RMSE = 13.9 m in Table 2.When conducting injection validation tests in the matrix, the model tends to underestimate the water level H at the monitoring points, with R 2 only 0.23953 and 0.22680, respectively.2f-h), similar differences are not observed-the flux fits well for the injection tests in the high-permeability fracture zone (Well2-2) and the matrix (Well3-1 or Well4-3).

K Estimates from the Head Data and Model Validation
Figure 3 shows the inversion results of Model_H based on the dataset of 15 water head monitoring points.In Figure 3a, the distribution of the high permeability fractured zone aligns well with the reference model in general, but the fractured zone's shape shows significant local differences.Particularly at x = 275 m, the fractured zone in the reference model narrows significantly, while in the inversion result, the fractured zone broadens instead.Additionally, although the inversion results depict the matrix on both sides of the fracture zone as low-permeability areas, the extent of these low-permeability areas is limited in the area where the observation wells are located, which does not match well with the reference model.This error corresponds to the uncertainty distribution shown in Figure 3b.Due to the lack of water head monitoring data, there are high-uncertainty regions around the upper/lower five monitoring points, indicating a significant error in depicting the low-permeability regions of the matrix.
As shown in Figure 3c-e for Model_H, based solely on the 15 water head monitoring datasets, water head H fits better for the injection test conducted at the monitoring point located in the high-permeability fractured zone (Well2-2) compared to the injection tests conducted in the matrix (injection points at Well3-1 or Well4-3), with R 2 = 0.38611, MAE = 9.28 m, and RMSE = 17.8 m in Table 3.When conducting injection validation tests in the matrix, the model tends to underestimate water level H at the monitoring points, with R 2 only 0.03772 and 0.04587, respectively.3f-h), similar differences are not observed-the flux fits well for the injection tests in the high-permeability fracture zone (Well2-2) and the matrix (Well3-1 or Well4-3).In Figure 4a, it can be seen that the addition of the flux data improves the depiction of the high-permeability fractured zone's shape and distribution, as well as the characterization of the matrix on both sides of the fractured zone.Especially for the characterization of the matrix on both sides, the flux information within the y = 100~400 m region reduces the uncertainty of the inversion results in this range (Figure 4b), making the distribution of low-permeability areas on both sides in Figure 4a   Figure 4c also shows results similar to Model_H, which used only the water heads in the inversion.The injection test conducted at the monitoring point located in the highpermeability fractured zone (Well2-2) fits water head H better with R 2 = 0.30778 than the injection tests conducted at the matrix monitoring points with R 2 = 0.14469 and 0.18428, respectively.There was no significant difference in the goodness of fit of the flux in the validation of the different injection tests.However, comparing the performance in Tables 3 and 4 with H-K and H-V-K shows that the addition of the flux data improves the validation of H, especially in the injection tests in the matrix zone (Well3-1 and Well4-3).

Differences between Head Data and Flux Data in K Estimation and Fractured Zones Characterization
This section will compare the inversion results using only water head data (Model_H, See Figure 3) and those using only the flux rate dataset (Model_V, See Figure 2), analyzing the impact of different data types on the inversion and on the identification and characterization of the high permeability fractured zone.
Figure 2a, which uses only the flux data in inversion, more accurately depicts the shape of the high-permeability fractured zone and also provides a more reasonable distribution of the low-permeability regions of the surrounding low-permeability bedrock.Additionally, the independent injection validation results for Model_V are generally better than for Model_H, especially in terms of the flux validation point (see Tables 2 and 3).
Although the inversion using flux data better characterizes the fractured zone distribution, the permeability value difference between the high-permeability region (K approximately 5 × 10 −3 m/s) and the low-permeability region (K approximately 5 × 10 −6 m/s) in the inversion results is only about three orders of magnitude, which fails to reflect the approximately five orders of magnitude permeability difference of those in the reference model.Although the estimated K is similar in the high-permeability zone, the estimated K of Model_H in the low-permeability region (~1.00 × 10 −7 m/s) is closer to that in the reference field; thus, it better reflects the significant permeability differences in the reference field.
Therefore, in the inversion of K to identify high-permeability fractured zones, the inversion model using only the flux data performs better overall in characterizing the fracture shape, but it is slightly worse in depicting specific K values.
It should be noted that in the above discussion, the number of flux monitoring points is greater than the number of water head monitoring points, which is probably the case in a real site.The increase in the number of monitoring points means that the informationgathering area is enlarged.The densely distributed flux monitoring points near the fractured zone enables the inversion model to reflect the distribution of the fractured zone more reasonably and accurately, thereby significantly enhancing the identification and characterization of the fractured zone.However, even if the number of flux monitoring points is several times that of the water head monitoring points, the inversion model's reliance solely on flux data does not necessarily guarantee more accurate permeability values compared to the inversion model's reliance solely on water head data.

Enhancements in Hydraulic Tomography through Flux Data Integration
The permeability distribution inverted from the data of 15 water head monitoring points (Model_H) generally reflects the position of the fractured zone (Figure 3a).Additionally, due to the limitation of the water head monitoring point, there are some errors in depicting the distribution of the low-permeability matrix.In Figure 5, when the number of water head monitoring points is reduced to five (Model_5H, see Supplementary Materials SI 1.1, Figure S1) or three (Model_3H, see Supplementary Materials SI 1.2, Figure S4), the inverted permeability distribution can still generally describe the position of the fractured zone, but there are significant errors in describing the whole.This indicates that when the input data is only water heads, the inverted permeability distribution is hard to depict in the fractured zone.This indicates that when the number of water head monitoring points is small, the inverted permeability distribution lacks sufficient information to constrain the iterative results, leading to non-uniqueness in the inversion results.The inversion results do not actually reflect the true permeability distribution.
The introduction of flux data has a significant positive impact on the HT inversion results.As shown in Figures 3a, 4a and 5, after introducing the vertically distributed flux data for joint inversion, the permeability distribution (Model_H-V, Model_5H-V, and Model_3H-V) is significantly improved, better depicting the position, the shape of the fractured zone, and the surrounding low-permeability bedrock.Even when compared to Model_H, which can roughly depict the position and distribution of the fractured zone, the introduction of flux data optimizes the inversion K distribution, the depiction of the fracture zone's shape, and the surrounding matrix distribution.Especially when the number of water head monitoring points is gradually reduced to five and three, the improvement in depicting the fracture shape with the flux data is more significant.
The above results indicate that introducing distributed flux data for joint inversion can effectively enhance the ability of the HT inversion method to identify and characterize the fractured zone.Previous research on the joint inversion of water head and flux data [20] shows that flux and water head data at the same monitoring point cannot be substituted for one another.Even with identical numbers and locations of monitoring points, HT inversion results using only flux data differ from those using only water head data.This indicates that flux data contains unique, non-redundant information.
This study confirms that the non-redundant information inherent in flux data plays an active role in the HT inversion process, particularly in characterizing the permeability structure.Because flux is vector data, it includes not only specific flux values but also flow directions.Velocity data is more sensitive to relative differences in permeability compared to specific permeability values, which enhances the depiction of permeability distribution.However, this characteristic also leads to poorer performance in estimating absolute permeability values.

Application Strategy of Head and Flux Data in Fractured Zone Characterization Using Hydraulic Tomography
Increasing the number of monitoring points is crucial for improving the results of the joint inversion of water head and flux data.The improvement due to increased monitoring points is more significant in joint inversion compared to using a single data type [20].However, deciding whether to increase head or flux data is a strategic challenge in the practical application of joint inversion to characterize fractured zones.According to the inversion model results, using data from 15 water head monitoring points can achieve a reasonably accurate permeability distribution (Model_H).However, this high spatial resolution, with five adjacent boreholes each containing three multilevel water head monitoring points, is costly and uncommon in actual field sites.
Therefore, considering the cost of the double-packer multi-level monitoring system for deep boreholes in bedrock, a scheme with fewer water head monitoring points, such as Model_5H or Model_3H, combined with a significant increase in vertically distributed flux monitoring points, would be more beneficial for practical field applications.Vertically distributed flux data can effectively constrain the non-uniqueness of the inversion model and accurately characterize the site's permeability structure.Additionally, the cost of obtaining flux data using the active DTS method is relatively low: the cost of one multi-level monitoring system in a single deep borehole could support the active DTS testing of all the wells at a site, making it advantageous for the practical application of the joint inversion method.
Considering that different types of monitoring data may vary in magnitude and dimensions, allocating weights to different types of observation data in joint inversion is also a challenge.Typically, the algorithm assigns equal weight to different types of data during the iteration process.When the magnitudes and dimensions of all observation data types are similar, this approach reflects the information of all data in a balanced manner.However, significant differences in magnitude or dimensions among different data types can result in the information contained in the smaller magnitude data being overshadowed by the larger magnitude data.Generally, the magnitude of flow flux data is often smaller than that of water head data.For instance, in this study's reference model, the magnitude of the flow flux data, about 10 −4 ~10 −9 , is much smaller than that of the water head data, which is generally greater than 10 2 .Directly incorporating both flux and head data into the inversion algorithm may lead to permeability estimation being predominantly controlled by the head data and thereby neglecting the flux data.
To solve this problem, the weight matrix in the algorithm can be optimized based on the differences in magnitude and dimensions among the observation datasets.However, considering that dataset differences vary by site, manually optimizing the weight matrix is costly.From a practical application perspective, this study proposes the following strategy for joint inversion methods of water head and flow flux data in fractured rock: First, perform the inversion based on dense flux data to obtain a permeability distribution with structural features as prior information.Then, combine this with limited water head data for joint inversion, optimizing the permeability values and distribution of the fractured zone.

Conclusions
This study investigates the potential of using HT joint inversion based on groundwater head and flux data to characterize fractured zones in bedrock, especially considering the impact of cost-induced differences in the quantity of head and flux data in actual field conditions.Synthetic injection tests were conducted in a fractured zone reference model to obtain water head and flux datasets of different quantities.Then, K distribution was estimated by hydraulic tomography, using the SimSLE based on pure water head data, pure flux data, and head-flux joint data.Additionally, validation tests were conducted in both the estimated permeability model and the reference model to obtain water head and flux at each monitoring point in the independent boreholes.By comparing the simulated water heads/flux and the uncertainty distribution of the estimated fields, the impact of different input data on the identification of fracture zone permeability was discussed.
By comparing inversion results with different sets of head data, flux data, and a combination of the two, it can be seen that the joint inversion of head and flux data can effectively enhance HT's accuracy in identifying the permeability characteristics of fractured zones.Especially when water head data is limited, which is often the case in actual conditions, flux information can significantly improve the identification and characterization of the permeability distribution.This is not only because the flux data itself brings non-redundant information but also because an increase in the number of monitoring points reduces the uncertainty of the inversion results.
Considering that the rapid development of active DTS method could provide highresolution spatial and temporal data on groundwater flux at relatively low cost, from the perspective of the practical application of HT in the characterization of fractured zones, it is better to estimate the permeability distribution base on the flux data, use it as prior information in further inversions, and combine them with limited water head data to optimize the permeability distribution.The findings of this study provide a promising method for determining the permeability distribution in real fractured rock sites and offer a relatively low-cost and efficient application strategy for detecting and characterizing fractured zones.

Figure 1 .
Figure 1.Vertical cross-section of the reference model and locations of head monitoring points and velocity monitoring points in the model; K = 5 × 10 −4 m/s is in the blue fracture area, and K = 5 × 10 −9 m/s is in the other area.Five monitoring boreholes (Well1-Well5) are placed within the X = 200 m to 300 m range, each spaced 25 m apart.All five boreholes penetrate the fractured zone.Due to the high cost of double-packer observation equipment, multiple water monitoring layers are limited in practice.Thus, each borehole has three water head monitoring locations at y = −100 m, −250 m, and −400 m, totaling 15 head monitoring locations.As mentioned in the the calculation to ensure the stability of the solving process.The covariance matrix r dd ε and the cross-covariance matrix r dy

Figure 2 .
Figure 2. Inversion results of Model_V and its validation.(a) is the inversion result based on the flux dataset; (b) is the uncertainty distribution with white circles indicating flux monitoring points; (c-h) are head and flux validation.To validate the performance of the K distribution obtained from the inversion, three new injection wells were activated in both the inverted K field and the reference model, as described in Section 2.1.2.The injection well locations and injection rate in the models are the same.The water head (H) and flow flux (Vx and Vy) distributions for each injection test

Figure 3 .
Figure 3. Inversion results for Model_H and its validation.(a) is the inversion result based on the head dataset; (b) is the uncertainty distribution with blue squares indicating water head monitoring points; (c-h) are validation scatter plots.

Figure
Figure 4a,b show the results of Model_V-H, combining the head data with the flux data.In Figure4a, it can be seen that the addition of the flux data improves the depiction of the high-permeability fractured zone's shape and distribution, as well as the characterization of the matrix on both sides of the fractured zone.Especially for the characterization of the matrix on both sides, the flux information within the y = 100~400 m region reduces the uncertainty of the inversion results in this range (Figure4b), making the distribution of low-permeability areas on both sides in Figure4amore reasonable.
Figure 4a,b show the results of Model_V-H, combining the head data with the flux data.In Figure4a, it can be seen that the addition of the flux data improves the depiction of the high-permeability fractured zone's shape and distribution, as well as the characterization of the matrix on both sides of the fractured zone.Especially for the characterization of the matrix on both sides, the flux information within the y = 100~400 m region reduces the uncertainty of the inversion results in this range (Figure4b), making the distribution of low-permeability areas on both sides in Figure4amore reasonable.

Figure 4 .
Figure 4. Inversion results of Model_H-V and its validation.(a) is the inversion result combining head and flux data; (b) shows the uncertainty distribution of the inversion results, with white circles indicating the input head monitoring points and black dots indicate the flux monitoring points; (c-h) are the validation scatter plots.

Figure 5 .
Figure 5.Comparison of inversion K fields using a different number of water head points.(a) is the inversion result based on 5 head monitoring points (Model_5H); (b) is the inversion result based on 5 head monitoring points and flux data (Model_5H_V); (c) is the inversion result based on 3 head monitoring points (Model_3H); (d) is the inversion result based on 3 head monitoring points and flux data (Model_3H_V).Blue squares indicate water head monitoring points and white circles indicate flux monitoring points.

Table 1 .
Details of synthetic injection tests in reference models with fault zone.
However, in the validation scatter plots for horizontal flux (Vx) and vertical flux (Vy) (Figure