Verification of Fracture Reorientation and Analysis of Influence Factors in Multiple Fracturing Treatment

A fracture will be initiated and propagated along the direction of maximum horizontal stress in fracturing treatment; however, in refracturing stimulation, the new fracture may be initiated and propagated along a different direction from the initial one. This is defined as a fracture reorientation. It is difficult to predict fracture reorientation due to the variation of formation properties after long-term production. To verify the existence of fracture reorientation and analyze its influencing factors in multiple fracturing treatment, experimental and numerical simulations are presented in this paper. Firstly, multiple fracturing stimulation is carried out with a self-assembled large true triaxial apparatus, and the fracture reorientation is successfully induced by changing the injection pressure and initial stresses in multiple fracturing processes. Then, numerical coupled hydromechanical modeling of the actual field production and injection well pattern is performed. In particular, the stress reversal region, which indicates the distance of fracture reorientation, and the factors that influence the reorientation are analyzed. The laboratory experiment and numerical simulation results show that the fracture reoriented angle obtained can be perpendicular to the initial fracture. Stress field and formation pressure are the two main factors that influence the fracture reorientation. With higher pressure differences and lower initial horizontal stress differences, the area in which it is possible to initiate reoriented fracture will be larger. The fractures of wells in the early production stage are hard to reorient due to the high formation and borehole pressure difference, and the fracture reorientation area will be expanded until the pressure difference is low to a certain value. This research result can guide oilfield stimulation treatments.


Introduction
Unconventional reservoir formations are considered to be with low permeability and porosity, and underground liquid seepage is hard to achieve during the producing process. Profitable productivity can be obtained only by implementing stimulation treatment [1,2]. However, after long-term production, the productivity of wells with initial fracturing treatment will decrease due to improper or small-scale fracturing treatment, closure of initial fracture, fines plugging, deformation and cracking of proppants, and borehole contamination [3,4]. In this situation, refracturing treatment is an effective stimulation method for wells to regain productivity. Initial fractures will be reopened, surface flushed, or proppant refilled in multiple fracturing treatment, and flow capacity will then be restored or even enhanced by extending the geometry of initial fractures. Additionally, by applying proper fracturing materials and large-scale fracturing treatment, fracture reorientation will be induced and new fractures will be initiated and extended to the direction different from the former one. This is the best scenario in field treatment because the fractures will penetrate regions with more residual oil and higher pressure.
Fracture reorientation is monitored in multiple fracturing stimulation in field treatments. Many researchers have studied the influencing factors that lead to fracture reorientation. Stress field changes are considered to be the leading factor of fracture reorientation in refracturing treatment [5]. By analyzing pressure data of refracturing field cases, it was shown that stress field changes due to propped initial fractures and depleted porous pressure is the main reason for fracture reorientation [6,7]. With theoretical analysis, it is concluded that the in-situ stress field would be altered in the production process from an existing initial fracture, resulting in new fractures that would be initiated and propagated to the direction perpendicular to the former one [8]. The initial reservoir stress field is steady before development; however, this status will be changed after well drilling, fracturing, liquid injection, and withdrawal. Regions of stress reversals caused by stress field alternation were observed in both vertical wells and horizontal wells [9,10]. For complicated well patterns, well injection and production make the stress field alternation difficult to predict. It was pointed out that fracture reorientation was always common in water flooding reservoirs due to the change of formation pressure after water injection [11]. It is difficult to evaluate stress field changes for those wells with complicated interference due to dynamic pore pressure change in the production and injection processes. With the development of fracturing technology, results of field and laboratory cases indicate that fracture reorientation can be induced by proper fracturing techniques and methods [12]. It has been verified that fracture reorientation can be easily induced by setting the direction of the perforation and controlling the horizontal stress difference in well completion and stimulation process. This can be done while the fracture initiation and propagation can be artificially guided to certain directions [13].
Mechanisms of fracture reorientation are identified and classified since it has become a research priority in restimulation treatments. For a fractured well, liquid production and injection will cause formation pressure change, and then, the stress field will be altered due to the increase or decrease of porous pressure. Meanwhile, it has been shown by numerical models that the stress field is an important factor that influences pressure and production performance [14]. Both the initial minimum and maximum horizontal stresses will be changed to some extent that the minimum can exceed the maximum one. The induced stress change will finally form an area where the directions of the initial stresses are totally reversed. This area is defined as the stress reversal region where fracture reorientation will exist. New fractures will be initiated at the direction perpendicular to the initial fractures in a reversal region [15]. Basic theories and models are applied to determine the mechanisms and influencing factors of fracture reorientation in multiple fracturing treatments.
In 1995, fracture reorientation caused by stress field change induced by production was first presented and explained based on the poroelasticity theory [16]. Since then, many models have been presented to deeply study the mechanisms and influence factors of fracture reorientation. Analytical and numerical models were built to calculate the stress changes from wells after fracturing treatment, and then, the influencing factors of the fracture length were discussed in refracturing treatment [17]. It was verified that there existed a stress reversal region around infilled wells after researching the stress changes induced by water flooding wells [18]. Coupled models were built and in situ stress changes due to well production were investigated, and models were applied to determine the right time for refracturing treatment. Meanwhile, factors that influence the existence of fracture reorientation were analyzed [19]. A heterogeneous, two-dimensional coupled model was presented to study the stress field changes due to production, and geological parameters that influence stress reversal were discussed [20]. Numerical models were presented to verify the possibility of fracture reorientation, and the relationship between enhanced oil recovery and fracture reorientation region was also discussed by Benedict and Wegner [21].
Based on a literature review, many theoretical studies and field tests have been done to study the mechanism of fracture reorientation in refracturing treatment. However, due to the difficulty of restoring reservoir conditions and simulating the complicated processes of multiple fracturing treatment, there are fewer researches on physical simulation in the laboratory to achieve fracture propagation. Therefore, fracture propagation in refracturing treatment can only be monitored by indirect ways in field tests, which are difficult to visually observe.
In this paper, the whole process of multiple fracturing treatment is completed on a self-assembled true triaxial apparatus that can apply three independent confining stresses, and fractured samples are sliced to directly observe the complicated fracture propagation. Then, numerical hydraulic-mechanic coupled modeling for a production and injection well pattern is simulated based on the poroelasticity and porous medial fluid theories. The main factors that impact fracture reorientation are further studied based on numerical models. Research results may guide field engineers in predicting the possibility of fracture reorientation before operating field treatment.

Fracture Reorientation in Laboratory Tests
Experiments were performed in the laboratory using a true triaxial apparatus, and multiple fracturing treatment was completed under three independent confining stresses applied to simulate realistic underground conditions. By applying various confining stresses and injecting fracturing liquids at different pressures, fracture reorientation is successfully induced and can be directly observed by slicing the test samples after finishing the experiments.

Test Samples and Equipment
2.1.1. Test Sample. The 30 × 30 × 30 cm 3 cubic samples were an artificial sandstone made from nonshrink grout/cement and sand that passes the # 40 or 0.42 mm opening mesh. The cement and sand were mixed with water in a concrete mixer and placed in a cubic mold to cure and dry. Two cubical samples were prepared and hydraulically fractured under 2 Geofluids different stages. The properties of the artificial sandstone are shown in Table 1, where ρ is sample density, ϕ is the porosity, k is the permeability, UCS is the unconfined compressive strength, BTS is the Brazilian Tensile Strength, and E is the Young's modulus. Tap water and hydrate pump oil were used as fracturing fluids.
2.1.2. Test Equipment. The true triaxial apparatus used in the laboratory simulation of fracture reorientation in cubical samples was designed and fabricated at Colorado School of Mines specifically for fracturing simulation studies. The apparatus is composed of five parts, namely, the true triaxial cell, hydraulic pumps to apply the three principal confining stresses, a digital data acquisition system, an acoustic emission monitoring system, and a miniature drilling rig. The main part is the true triaxial cell which consists of a rigid cylindrical frame with rigid plate lids on top and bottom. The confining pressure is applied via three dog-bone flat jacks hydraulically pressurized by manual pumps and are resisted by three rigid plates at the opposing sides of the flat jacks. Fracturing is induced through a miniaturized casing to a bottom-hole well that is partly uncompleted. Fracturing pressure is applied by two Isco micropumps working in tandem. Figure 1 shows a schematic diagram of the assembled true triaxial fracturing system with multiple well locations. Detailed information on the design and operation of the true triaxial system can be found in Lu et al. [22] and Frash et al. [23].

Experimental
Procedures. The experimental procedures to simulate fracture reorientation in rocks include four main steps for one simulation stage: (1) sample preparation and equipment assembly, (2) simulation of reservoir conditions by applying confining pressures to the three sample faces, (3) simulation of the fracturing treatment by increasing pressure in the injection well using the Isco pumps, and (4) posttest analysis by slicing the samples into thin cuts and imaging the induced fracture(s). The fracture initiation and propagation processes are monitored in real-time from the injection pump pressure and by acoustic emission sensors. Initial stresses applied to the samples and locations of the acoustic sensors are shown in Figure 2. S 1~S6 are the six AE sensors; σ h , σ H , and σ v are the minimum horizontal stress, maximum horizontal stress, and vertical stress, respectively. Five boreholes were drilled in sample 1 under the following test procedures: (1) the vertical, maximum horizontal, and minimum horizontal stresses were set to 16, 12, and 8 MPa, respectively. Then, the middle well was fractured with water at a flow rate of 0.5 mL/min. The injection pump is stopped when a pressure drop was observed, which indicated that fracture has been initiated, or when the acoustic emissions indicated fracture initiation. (2) The flow rate was increased to 30 mL/min to refracture the sample. (3) The fracturing fluid is injected into the four edge wells with a constant bottom-hole pressure to prevent the fracture from extending to the corner wells. A high-viscosity oil is then injected into the middle well, and the pump was kept running until the fracture propagated close to the edge of the sample. (4) The major and minor horizontal stresses were rotated 90°, and the middle well injection flow rate was increased to 20 mL/min. Figure 3 shows the pressure change with time during the three-stage fracturing process.
Three boreholes were drilled in sample 2 using the following procedures: (1) the three principal stresses were set to the same magnitudes as in sample 1. Then, the corner wells were pressurized at a constant flow rate of 0.5 mL/min. The pump was kept running until fractures were fully propagated.
(2) The injection fluid was changed to a high viscosity fluid. The middle well was fractured with a high constant flow rate of 30 mL/min, and the pump was stopped as soon as fracture initiation occurred based on a well pressure drop or acoustic emissions; then, the well pressure was released. (3) The fracturing fluid was injected into the middle well at constant pressure to maintain a constant bottom-hole pressure, and then, another corner well is fractured. Figure 4 shows the pressure change with time during the multiple fracturing process.

Experimental Results and Analysis.
After the fracturing tests, the samples were sliced into thin horizontal sections about 2.75 cm in thickness to observe the induced fractures directly. The sliced sections showed direct evidences of fracture reorientation induced by changing the initial stress field and the injection pressure. For sample 1, the fracture initiated in the first fracturing treatment and then propagated in the secondary fracturing treatment. By rotating the orientations of the initial horizontal stresses, a second fracture initiated at the direction perpendicular to the initial fracture and then propagated as shown in Figure 5. Fracture propagation was also located by acoustic emission as shown in the 3D scatter diagram, which shows a high correspondence of the AE signals with the orthogonal-induced fractures. Spots of different colors show the fracture location of every fracturing treatment. Backpressure applied to the corner wells showed that the fracture intersected the corner wells.
The acoustic emission events for fracture location in the sample during multiple fracturing indicate fracture initiation and propagation of every fracturing treatment, which agrees with the observations from the sliced samples.
As shown in Figure 6, multiple tiny fractures are initiated around the middle well, and the initiated directions of these fractures are complicated by fracturing at a high flow rate. To observe the fractures, pressure is released as soon as possible after fracture initiation to avoid the main fracture intersecting with the tiny fractures. Fracture propagation located by acoustic emission is shown by the 3D scatter diagram. With the same stress field, well 2 was fractured before well 1, and the fracture propagates along the initial maximum horizontal stress. However, the fracture propagation of well 1 is different with that of well 2 because of the applied 3 Geofluids backpressure of well 1 and altered stress field. It is a complicated process because the fracture initiated and the propagated direction can be impacted by so many factors.

Fracture Reorientation Mechanism
Analysis. It can be verified from the experiments that fracture reorientation occurred in the process of refracturing operation, and new fractures could initiate at angles different from that of the initial one. In refracturing treatment, injected pressure, bottomhole stress distribution, and well interference are all important factors that determine the initiated direction of new fractures. In the refracturing process, the pore pressure around the bottom well will be increased sharply with highpressure liquid injected. However, due to the existence of an initial fracture, the pore pressure at the parallel direction with the initial fracture is increased more than that at perpendicular direction, which leads to an unbalanced increase of the minimum and maximum horizontal stresses. As shown in Figure 7, the inner stress reversal region is formed which determines the initiation direction of fracture reorientation.
However, it is not enough for field applications only to have new fractures initiated and reoriented. It is important that fractures must propagate and extend to the far field away from the well bottom. The propagation distance of fracture reorientation is determined by the stress fields around and far away from the well bottom. This stress reversal region is determined by the stress field far away from the well bottom-hole which is induced by liquid production and injection, and the outer elliptical stress reversal region determines the extended distance of fracture reorientation as shown in Figure 7.
Based on the above interpretation of the experimental results and the fracture reorientation mechanism, it can be shown that new fractures will be initiated around the well bottom-hole during the process of refracturing treatment. By applying a high flow rate of injected liquid and changing the in situ stress field in a small-scale sample in laboratory tests, it is easier to create new reoriented fractures. Oriented perforation is one common technique in field treatment, but it is difficult to simulate in small-scale cubic samples in the laboratory, especially in refracturing treatment. However, it can be inferred that with oriented perforating technology in field treatment, fractures will be easily reoriented in the initiation period. As analyzed above, if the outer stress reversal region is small or even nonexistent, new fractures will propagate along the original direction. In order to identify the distance of the reoriented fractures, it is necessary to study the far stress field distribution which will be varied due to longterm production and injection, which is known as stress variation induced by a poroelastic effect. Since it is impossible to simulate samples that are large enough for poroelastic effect research, numerical simulation is used in the following to extend the experimental results.

Coupled Hydromechanical Modeling
As mentioned above, the formation of liquid withdrawal and injection will lead to pressure change and will cause variation of the stress field; the formation stress field determines the direction of fracture initiation and propagation in refracturing treatment. Meanwhile, the area of the reversal stress region determines the distance of fracture propagation. Therefore, research on the stress field induced by pore pressure is an important task before multiple fracturing treatments. To furtherly study the induced stress distribution after initial fracturing treatment, numerical coupled hydromechanical models of production and injection well patterns are simulated to study the factors that influence the stress reversal region.
3.1. Poroelastic Model. The deformation of porous media saturated with liquid is controlled by the applied effective stresses. The relationship between the strain tensor and the volumetric strain for porous media can be expressed as:  Geofluids where ε is strain tensor of porous media, u is the displacement of porous media, ε v is the volumetric strain for porous media, and trðεÞ is the trace of the strain tensor ε. For linear poroelastic media, based on the effective stress principle, the relationship of the effective strain, the overall stress, and the pore pressure can be expressed as: where σ ′ is the effective stress tensor (in MPa), σ is total stress tensor (in MPa), p is the pore pressure of the porous media (in MPa), a is Biot's poroelastic constant, and I is second-order identity matrix. The momentum equation for the motion of elastic porous media can be obtained from equilibrium: When the rock reaches the state of mechanical equilibrium, the velocity of the point of mass is zero, so the first expression of the equation is zero. By combining Eq. (3), the above equation can be expressed as: Then, associating Eq. (1) and Eq. (2) with the above equation yields: where G is the shear modulus of porous media, G = E/ð2ð1 + vÞÞ (in MPa), v is Poisson's ratio, and λ is the Lamé constant (in MPa) defined as λ = Ev/ð1 + vÞð1 − 2vÞ.

Model for Fluid Flow in Porous Media.
In this paper, the hydraulic response of a reservoir is modeled as a two-phase immiscible fluid flow. The model is applied to reveal pore pressure changes in the reservoir during fluid injection and production processes. Considering the deformation of the porous media during water flooding, the kinematic velocity of the porous media skeleton and the actual velocity of the pore fluids in porous media is the sum of the flow velocity of liquid and deformation velocity of porous media.
where U is the fluid seepage velocity in the porous medium (in m/s), U r is the relative velocity between liquid and porous medium (in m/s), v s is the velocity caused by    6 Geofluids the deformation of the porous medium (in m/s), v is the Darcy seepage velocity of the fluid (in m/s), and s j is the saturation of the liquid phase j. Therefore, the absolute velocity of liquid phase j in the coupled hydromechanical models can be expressed as: Based on the law of mass conservation, the equation of continuity for the two-phase oil/water seepage can be expressed as: where B j is the volume factor of liquid phase j (in m 3 /m 3 ) and Q j is the flow rate of the injected or produced liquid at standard condition (in m 3 /d).
Combining Eq. (8) and Eq. (9) yields: where μ j is the viscosity of liquid phase j (in mPa·s) and λ j is the mobility of liquid phase j (in 10 -3 μm 2 /(mPa·s)) and is defined as λ j = kk j /μ j . The above equation can be expanded as: The compressibility factors of the rock pore volume, oil, and water and can be defined, respectively, as:

Geofluids
From the continuity equation for the porous medium, the relationship between the strain and displacement can be obtained as: Combining Eq. (11), Eq. (12), and Eq. (13) with Eq. (10) and multiplying by B j give the saturation equation of liquid phase j: The pressure equation can be obtained by combining the oil phase of Eq. (10).
The relationship of oil saturation and water saturation can be shown as: The coupled hydromechanical model described by the above equations can be discretely solved by a finite element method.

Coupling of Stress and Seepage in the Artificial Fracture
Region. The relationship between the fracture width and stress can be expressed as [24]: where w f is the width of the artificial fracture (in m), w fi is the initial width of artificial fracture (in m), K n is the normal stiffness of artificial fracture (in MPa/m), and σ n ′ is effective normal stress acting across the fracture surface.
Based on the seepage model of two flat plates, the relationship between the fracture permeability and fracture width can be expressed by using the parallel-plate law as: A dimensionless fracture conductivity R conductivity can be defined as the ratio of real-time fracture conductivity and initial fracture conductivity.

Coupling of Stress and Seepage Field in the Rock Matrix
Region. The porosity of the porous medium is defined as: where V is the total volume of porous media (in m 3 ) and V s is the skeleton volume of porous media (in m 3 ). Combining with volumetric strain, dynamic porosity can be defined as: Based on the capillary bundle model, dynamic permeability can be defined as: 3.3.3. Boundary Conditions. The boundary conditions of the poroelastic model include the displacement boundary conditions and stress boundary conditions. The displacement boundary conditions at the initial time of the reservoir can be shown as: The stress boundary conditions at the initial time of the reservoir can be shown as: where n is the stress in the normal direction,

MPa.
Based on the seepage theory, the boundary conditions for the seepage model include the initial time conditions and boundary conditions. The initial time conditions of the seepage model can be shown as: The interior and exterior boundary conditions include the Dirichlet and Neumann conditions. The Dirichlet boundary condition is applied when the liquid pressure is known at the boundary, and the Neumann boundary condition is applied as the liquid pressure gradient at the normal direction of the exterior boundary. neglecting the capillary force and the gravity, neglecting the temperature change, and only linear elastic deformation for the rocks and without the fracturing occurring during the producing process. COMSOL Multiphysics is applied for discretizing and implementing the coupled hydromechanical model [20]. All data used in the simulation were collected from a typical low permeability reservoir, and the basic parameters used in the model are shown in Table 2.
The rhombus inverted nine-spot well pattern is widely applied in the development of low permeability reservoir, and for the convenience of numerical simulation and research on the well interference of target wells, this approach is played in choosing wells for numerical simulation, as shown in Figure 8. Two wells with different locations are set as target wells for detailed research; P3 is an edge well and P5 is a corner well.

Model Results.
To investigate the formation stress variation of fractured vertical wells in the water flooding reservoir, the models of the well pattern are simulated for producing 3000 days. Change of maximum horizontal stress, in situ stress, and strains at the direction of the initial minimum and maximum horizontal stresses (x and y axis) are analyzed with different field development periods. Figure 9 shows the maximum stress distribution in different development periods. Figure 10 shows the minimum stress distribution in different development periods. Figure 11 shows the shear stress distribution in different development periods.
By analyzing the stress distributions in different locations, it is found that stress changes are mainly around oil wells and water wells. The low-stress region is formed around fractured oil wells due to the pressure drop caused by liquid production, while the high-stress region is formed around water wells due to high pressure caused by water injection. The stress gradient is high around fractured oil wells and water wells. By analyzing the stress distribution in different development periods, stresses change with development, which increases around water wells and decreases around oil wells.
The maximum stress directions in different development periods are shown in Figure 12, and the direction at the initial time is along the propagation of the fracture. In the early period, the direction of stress varies sharply around all oil wells, the changed direction of stress is reversed to be perpendicular to the initial direction, and the stress direction is changed to the radial direction of the wellbore. In the midperiod, the stress reversal region expands while the area of the middle well is bigger than that of the edge wells. In the late period, the stress change around water wells is dominant and expands at the radial direction of the wellbore, the stress reversal region decreases and then gets steady for the edge wells, while it increases and then gets steady for the middle well.

Parametric
Analysis. The effects of the initial in situ stresses and fluid pressures, which are considered to be the dominant factors that influence the fracture propagation and stress reversal region, will be analyzed separately in this paper.

Initial Horizontal Stresses.
To investigate the influence of stresses on fracture reorientation, the models with different ratios of maximum and minimum horizontal stress of 1.02, 1.04, 1.06, and 1.12 are simulated. In this paper, the dimensionless fracture reorientation factor (L a ) is defined as the ratio of the equivalent radius of the reversal region and half-length of the initial fracture. The results are shown in the following figures. Figures 13 and 14 show dimensionless fracture reorientation factor changing with different initial horizontal stress ratio for the edge well and the middle well, respectively.
It can be seen that for edge wells, the dimensionless fracture reorientation factor increases to a maximum and then decreases with time, while decreases with the increasing initial horizontal stress ratio. It means the initial horizontal stress ratio has less influences on the time when the stress reversal region gets to be steady, while influences more on stress reversal distance. For the middle well, the dimensionless fracture reorientation factor increases and then slows down with time, and it also decreases with the increasing initial horizontal stress ratio. A small horizontal stress ratio is more conducive to fracture reorientation.   Figures 15 and 16 show the dimensionless fracture reorientation factor changing with different production pressure difference for the edge well and middle well, respectively.
For the middle well, with the same bottom pressure, dimensionless reversal distance increases with higher injec-tion pressure, when the bottom pressure is lower than 25 MPa. The reservoir pressure depletes, and the pressure gradient is small in the development process, which causes that dimensionless reversal distance to increase and then decrease gradually. For the edge wells, dimensionless reversal     distance increases and then decreases with time. The same parameter increases with higher injection pressure, which means with higher producing pressure difference, stress variation induced by production is bigger, and fracture reorientation is easier due to high horizontal stress difference.
With the same injection pressure, for the middle well, dimensionless stress reversal distance increases with high producing bottom pressure at the same development period. With higher bottom pressure over 8 MPa, it is difficult for stress reversal to occur due to the low pressure difference around the producing well (with initial reservoir pressure of 9 MPa), and the stress reversal region expands until the production time of 900 days. For the edge well, the dimensionless stress reversal distance increases and then decreases with time. The same parameter increases with the decreasing bottom pressure of the producing well. With the higher bottom pressure of over 8 MPa, it is difficult for stress reversal until producing time of 100 days. The results indicate that the stress change for the edge well is faster than that of the middle well.

Conclusions
The mechanisms of fracture reorientation in refracturing treatment are studied, and the concept of the two stress reversal regions is proposed. With physical laboratory tests and numerical simulation, the fracture reorientation is verified to be existed around the well bottom, and influence factors that impacted fracture extended distance are analyzed.
(1) Two stress reversal regions that impacted fracturing treatment in refracturing treatment are proposed: the inner one is highly influenced by fracturing parameters during the refracturing process, and the outer one is impacted by the development of the reservoir

Geofluids
(2) Fracture reorientation is induced by changing the initial stress field and the injection pressure in laboratory tests on the true triaxial fracturing system. The direction of reoriented fracture can reach to be perpendicular to the initial fracture (3) In the refracturing treatment, new fractures are easily initiated and reoriented around the well bottom by controlling fracturing parameters, which is verified by experimental results (4) The propagation distance of the reoriented fracture is determined by poroelastic effect, the initial stress field will change with pore pressure varies, and the outer stress reversal region is formed which determines the fracture reoriented distance in the fracturing treatment (5) The stress field and the injected and produced pressures are the main factors that influence the stress reversal region, and the propagation of the reoriented fractures will be varied due to different well location

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.   12 Geofluids