An Analytical Solution for Variable Viscosity Flow in Fractured Media: Development and Comparative Analysis With Numerical Simulations

Explicit fracture models often use analytical solutions for predicting flow in fractured media, usually assuming uniform fluid viscosity for simplicity. This assumption, however, can be inaccurate as fluid viscosity varies due to factors like composition, temperature, and dissolved substances. Our study, recognizing these discrepancies, abandons this uniform viscosity assumption for a more realistic model of variable viscosity flow, focusing on viscous displacement scenarios. This includes instances like injecting viscous surfactants for hydrocarbon recovery in fractured reservoirs or soil decontamination. This presents a significant challenge, enhancing our understanding of transport within fractures, mainly governed by advection. Our study centers on a low‐permeability rock matrix intersected by two fractures with variable apertures. We employ two methods: an analytical approach with a new solution and numerical simulations with two distinct in‐house codes, discretizing both the rock matrix and fractures with two‐dimensional triangular elements. The first code uses a Discontinuous Galerkin finite element method, while the second utilizes a finite‐volume method, allowing a comprehensive comparison of solutions. Additionally, we investigate parameter identifiability, like fracture apertures and viscosity ratios, using breakthrough curves from our analytical solution, applying the Markov Chain Monte Carlo technique.


Introduction
Variable viscosity flow (VVF) arises from the interaction of fluids with different viscosities within a porous medium, resulting in intricate flow patterns (Chou et al., 2015;Yih, 1961).These variations in viscosity can arise from compositional differences, temperature gradients, or the presence of dissolved substances (Animasaun, 2015).A comprehensive understanding of VVF in porous media is of utmost importance in numerous applications, including groundwater contamination analysis, enhanced oil recovery, geothermal system simulations, and subsurface storage (Kim, 2014;Prasad et al., 2011;Strack, 2019;Zidane, 2023).
Numerical models developed to simulate VVF in porous media incorporate various mathematical formulations, including the Navier-Stokes equations, Darcy's law, and the continuity equation.These formulations often include additional terms to account for the spatial variation of fluid viscosity in the Navier-Stokes equations and Darcy's law (Badea et al., 2010), or through the density term in the continuity equation (Jayanthi & Kumari, 2007).Numerical modeling techniques, such as finite difference (Sabet et al., 2018), finite element (Caucao et al., 2017;Sabet et al., 2018), and finite volume (Schmid et al., 2013) methods, have been extensively utilized to solve the governing equations of VVF.Through numerical simulations, researchers have investigated the influence of viscosity variations on flow patterns (Reddy & Reddy, 2013;Shankar et al., 2017), pressure distribution (Zaytoon & Hamdan, 2021), and velocity profiles (Agrawal et al., 2020) within porous media.Furthermore, studies have examined the effects of various parameters, such as temperature, chemical reactions, and composition changes, on the behavior of variable viscosity fluids in porous media (Ellahi & Afzal, 2009;Graf & Boufadel, 2011;Kandasamy and Hashim, 2008).Analytical models have also been developed to provide insights into the fundamental aspects of VVF in porous media (Ananthaswamy & Uma Maheswari, 2015;Angot, 2021;Ellahi & Afzal, 2009;Hayek, 2019).These analytical solutions often focus on simplified scenarios and idealized geometries to derive mathematical expressions that describe the flow behavior.These models can offer analytical relationships between variables such as pressure, velocity, and viscosity (Al-Khafajy, 2016;Choudhary et al., 2023;Ellahi & Afzal, 2009;Suk, 2017), enabling researchers to gain a deeper understanding of the underlying mechanisms governing VVF.Furthermore, these analytical approaches are valuable in validating and benchmarking numerical simulations (Vanderborght et al., 2005), and provide approximate solutions in cases where complex numerical simulations are not feasible or necessary (Baigereyev et al., 2021).
The existing body of research primarily revolves around numerical modeling approaches to study VVF in fractured media.Although certain analytical models for flow and transport in fractured media have been proposed (such as the Generalized Integral Transform Technique (GITT), discrete fracture networks (DFNs), and stochastic continuum (SC)) (Ahmed et al., 2019;Cotta et al., 2020), it is noteworthy that the development of analytical solutions specifically targeting VVF in fractured media has received very limited attention.Therefore, the objective of this study is to develop a novel analytical solution that considers VVF in porous media with explicitly defined fractures.In doing so, our goal is to address the existing research gap and contribute to advancing the understanding of VVF.In addition, we aim to establish a benchmark solution that can be utilized for validating numerical models while providing efficient and approximate solutions to expedite VVF simulations.
The remainder of this paper is organized as follows.In the subsequent section, a review of related work on the development of analytical models for flow and transport in fractured media is presented.To establish a solid theoretical framework, Section 3 includes a review of the governing equations that describe VVF in fractures.The test case, assumptions, and development of the benchmark analytical solution are explained in Section 4. The numerical model employed to simulate the same test case is elaborated in Section 5. Section 6 presents and discusses the results obtained from a comparison between analytical and numerical simulations.Section 7 concludes the study and reviews the main findings.

Related Work
Flow and transport modeling in fractured media can be approached using either an implicit or explicit representation of fractures (Berre et al., 2019).In the implicit approach, fractures are accounted for using single or multi-continuum models.In a single continuum model, the equivalent permeability of a porous medium considers the properties of the fractures and fracture network (Durlofsky, 1991;Liu et al., 2016).By contrast, multicontinuum models utilize two or more superimposed media with distinct flow and/or transport equations (Fahs et al., 2014;Jourde et al., 2002;Kordilla et al., 2012).Alternatively, the explicit approach adopts an explicit representation of fractures in the modeling process.This approach directly represents individual fractures and accounts for their geometrical and hydraulic properties in the simulation.Explicit representation allows for a detailed characterization of the fracture network and captures the impact of individual fractures on flow and transport behavior (Hirthe & Graf, 2015).The incorporation of explicit fracture representation fosters more accurate modeling, encapsulating the behavior of fluid flow affected by both the characteristics of fractures and variations in fluid viscosity (Huang et al., 2022).
Examples of analytical flow and transport models that rely on the explicit representation of fractures include Tang et al. (1981), Wu et al. (2004), Xing et al. (2017), O'Malley et al. (2018), Cotta et al. (2020), andHyman et al. (2022).The mathematical framework for modeling flow and transport in fractured media is generally provided by a combination of equations from fluid dynamics (e.g., Darcy's law and the continuity equation), mass transport (i.e., advection-dispersion equation), and if relevant, heat transfer, as well as equations describing the specific properties of the porous media and fractures (Gao & Ghassemi, 2020;Hyman et al., 2022;Khoei et al., 2023).The latter describes the hydraulic and mechanical properties of fractures, such as their aperture, permeability and stiffness, and may involve empirical or theoretical relations (Wang et al., 2019).
When creating an analytical model that includes fractures, several simplifying assumptions are commonly adopted to make the model more manageable, albeit at the risk of introducing certain inaccuracies or errors.These Water Resources Research 10.1029/2023WR035741 assumptions are: (a) fractures are modeled as flat, reduced-dimensional entities, disregarding any curvature or unevenness in three dimensions (Hyman et al., 2022); (b) it's assumed that each fracture maintains a consistent width and transmissivity (Hirthe & Graf, 2015); (c) the model overlooks the potential impacts of diffusion through the surrounding rock or adsorption processes; and (d) The flow of fluid within each fracture is considered to be smooth and evenly distributed (Fang & Zhu, 2018).It's important to point out that in the analytical solution for a single fracture proposed by Tang et al. (1981), which was later utilized by Graf and Simmons (2009), matrix diffusion was indeed considered.Previous analytical solutions typically concentrated on specific fracture setups.For instance, Zhang and Ayala (2019) developed solutions for a single fracture with changing viscosity, while solutions for flow through parallel fractures were addressed by Zhou andWang (2023), andToller (2022) focused on intersecting fractures.
The uniform fluid viscosity assumption is commonly applied in modeling flow through porous and fractured media, as it simplifies the mathematical formulation and solution of the problem (Baigereyev et al., 2021;Berre et al., 2019).However, fluid viscosity can actually vary due to different factors (Reddy & Reddy, 2013).VVF is significant in some applications, such as in the injection of viscous surfactants for enhanced hydrocarbon recovery in fractured reservoirs or for soil decontamination.Neglecting viscosity variations can lead to significant errors in model predictions (Hasona et al., 2018).To enhance the accuracy of our models and to better represent real-world flow and transport in fractured media, this study aims to move beyond the usual assumption of uniform viscosity.This shift increases the complexity of our mathematical model as fluid properties now vary with space and time, representing a substantial departure from traditional models where viscosity is considered a constant.

Theoretical Background
In this section, the equations governing VVF in fractures are discussed.For simplicity, the flow and transport are considered for a horizontal fracture network.The flow of fluids with variable viscosity in fracture networks can be described by the steady-state mass conservation equation and Darcy's law, respectively (Simacek & Advani, 2003): where P is the pressure [Pa], q is the Darcy velocity [LT 1 ], μ is the fluid dynamic viscosity [ML 1 T 1 ], k is the permeability of the fracture [L 2 ], and ∇ represents the gradient operator.The mixing process of dispersion and diffusion is often disregarded within fractures, based on several assumptions: first, the relatively short length of the fractures; second, the absence of soil filling these fractures; and third, the relatively high velocity of fluid flow inside these fractures.Hence, the transport of contaminants can be assumed to be governed only by advection as follows: where V = q/θ is the fluid velocity [LT 1 ], θ is the porosity [-] and C is the nondimensional relative concentration [-].Equation 3 includes a time-dependent term to capture the transient dynamics of contaminant transport, enabling the analysis of scenarios where the contaminant is introduced into the system at a specific point in time.
The coupling between flow and transport equations occurs through density and viscosity relationships that are dependent on concentration (Ackerer et al., 1999).Considering fractures as open empty channels without a porous medium (θ = 1), the permeability of the fractures can be approximated based on the aperture of the fracture (e) [L] (Adler et al., 2013): which is often referred to as the cubic law (Witherspoon et al., 1980).

Developing the Analytical Solution
This section offers a summary of the case under examination, detailing the essential parameters, boundary conditions, and other pertinent information.Following sections delve into analytical solutions for flow modeling in this test scenario, considering both uniform and variable viscosities of the fluids being injected and displaced.

Description of the Test Case
We examine a rock matrix intersected by two saturated fractures, each of length L, but with differing thicknesses, which are joined at their ends, as depicted in Figure 1a.A viscous fluid is injected at the upstream node at a fixed flow rate.The initial fluid within the fractures is displaced by the injected fluid, exiting through the downstream node.This scenario is converted into a synthetic test case, illustrated in Figure 1b, featuring a domain with a rectangular shape measuring 10 cm in length and 50 cm in width.The domain consists of a rock matrix characterized by an extremely low permeability k m which is considered inactive in the context of fluid dynamics.This means that the flow of fluids through this part of the domain is negligible due to the minimal permeability characteristics of the matrix.
Within this domain, there are two parallel fractures of different apertures that cross the central zone from top to bottom.The first fracture is located at x = 4 cm and has an aperture of e 1 , whereas the second fracture, situated at x = 6 cm, has a smaller aperture of e 2 .The domain is initially free of contaminant.The boundaries of the domain are impermeable, except for the injection point at the top, located at x = 5 cm and y = 50 cm, where continuous injection introduces a fluid with a normalized concentration C inj = 1 at a specified inlet flux Q.The injected and displaced fluids exit the system through the outlet orifice at the bottom, positioned at x = 5 cm and y = 0.At the outlet, a Dirichlet boundary condition with zero pressure is applied.The injected fluid has a viscosity of μ 1 , whereas the displaced fluid has a viscosity of μ 0 .It is assumed that the two fluids do not mix inside the fractures (dispersion is neglected) and mixing between fluids occurs only in the output node.Notably, the viscosity contrast can be greater than the permeability contrast if the injected fluid is significantly more viscous than the displaced fluid.
During subsequent analysis, an analytical solution is formulated for flow and transport scenarios.Initially, the analysis focuses on a tracer case, where fluids with identical viscosities are injected and displaced.Subsequently, the analysis explores the scenario of viscous displacement, where the injected fluid possesses significantly higher viscosity than the displaced fluid.In terms of transport within fractures, the influence is attributed to advection processes, whereas the impact of diffusion processes within fractures is disregarded.

Analytical Solution for Uniform Viscosity of Injected and Displaced Fluids
Initially, an analytical solution is developed for the scenario where the fluids that are injected and displaced have the same viscosity.Noting that the total water flux Q [L 2 T 1 ] is the sum of the flux through the first and second fractures (denoted by Q 1 and Q 2 respectively), the mass conservation equation for the fluid can be expressed as follows: where V i = Q i /e i is the fluid velocity [LT 1 ] inside the fracture i = 1, 2. The mass conservation Equation 1, written for the fracture i with a constant aperture e i and filled by a fluid of constant viscosity μ 0 , implies that the pressure has a linear variation along each fracture.Assuming the pressure difference between the entry and outlet ΔP=P ent -P out [Pa], Equation 2gives: Using Equation 5and assuming that δ = e 2 /e 1 , we obtain: The fluid velocities V 1 and V 2 remain constant over time.The pressure remains steady over time and exhibits a linear variation along the x-axis, starting from P ent at x = 0 and transitioning to P out at x = L. Assuming P out = 0, the input pressure P ent can be calculated using Equation 6as follows: The advection Equation 3 implies that the concentration C is constant along the characteristic x(t) = X t; x, t) passing through the point x,t) and determined by the following system of equations: Which implies that C (x,t) = C (0,t*) where (x,t) backtracks to (0,t*).Thus, for a continuous injection with a constant input concentration C inj in a domain initially free of contaminant, the concentration in the fracture is either The mass conservation of the contaminant at the outlet reflects the ideal mixing of the fluids arriving from F1 and F2: where C L is the outflow concentration, C 1 (L,t) and C 2 (L,t) denote the concentrations at the end of the first and second fractures, respectively, each with a length of L. Hence, substituting Equation 7 into Equation 10 gives the outflow concentration: The breakthrough curve (BTC) illustrates the arrival of a rapid front from fracture 1 (F1) at the outlet at t 1 = e 1 (1 + δ 3 ) (L/Q) followed by a slower front from fracture 2 (F2) at t 2 = e 2 (1 + δ 3 ) (L/Q).Before t 1 , the outlet concentration is zero, while it becomes one after t 2 .Between t 1 and t 2 , the fluid exiting F1 has a concentration of 1, while the fluid exiting F2 has a concentration of zero.The resulting outlet concentration is determined by the mixing of these two fluids, as follows:

Analytical Solution for Different Viscosity of Injected and Displaced Fluids
In situations where the viscosities of the injected and displaced fluids are different, the flow rates within the fractures change over time, resulting in three distinct phases during the infiltration process.In the first phase, the leading edge of the injected fluid hasn't reached the outlet in either of the fractures, a milestone that's expected at time t 1 .In the second phase, the injected fluid in F1 arrives at the outlet, but it hasn't yet reached the outlet in F2, an event anticipated at time t 2 .Lastly, in the final phase, both fractures are fully saturated with the contaminant.The following sections will detail the analytical solutions for each of these phases separately.

The Initial Flow Regime (t = 0)
At the onset, both fractures are occupied by the initial fluid, which is free of contaminants and has a viscosity denoted as μ 0 .The initial velocities in this scenario (V 1 (0) and V 2 (0)), similar to those in the tracer case, are determined by Equation 7.

The Flow Regime During the First Phase (0 < t ≤ t 1 )
In the first phase, fracture i, which is L in length, contains the injected fluid with a viscosity of μ 1 in its segment X i (t) and the displaced fluid, having a viscosity of μ 0 , in the remaining segment (L-X i (t)).This configuration is illustrated in Figure 2.
We define P * i as the pressure at the boundary where the injected and displaced fluids meet within fracture i.Consequently, integrating conservation Equation 1 with Darcy's Equation 2 reveals that the pressure profile is linear in the first segment (0 ≤ x ≤ X i (t)) of the fracture, which is filled with the injected fluid, and similarly linear in the latter segment (X i (t) ≤ x ≤ L), occupied by the displaced fluid.The continuity of flux across these two segments within fracture i is expressed as follows: and Detailed calculations presented in Appendix A yield the position of the interface in F1, as detailed in Equation A14: Therefore, the fluid velocity in F1 is: The interface position in F2 is given by combining Equation 14and Equation A5: X The fluid velocity in F2 is obtained from Equation 5: The time t 1 when the contaminant in F 1 reaches the outlet (X 1 (t 1 ) = L) is obtained by solving the following equation: At this point in time, the position of the front in the two vertical fractures X 1 (t) and X 2 (t) are respectively given by Equation 14and Equation 16and the velocities V 1 (t) and V 2 (t) inside the fractures are respectively obtained from Equation 15 and Equation 17.
Denoting X 2 (t 1 ) = X 2 * , the position of the contaminant front in F2 at the end of the first phase (t = t 1 ) (i.e., when the contaminant in F 1 reaches the outlet).Using Equations 16 and 18, we obtain: At t = t 1 , the velocity in the two fractures is: Contrary to the tracer test case, the pressure does not remain constant over time.To determine the pressure distribution inside F1 (or F2) at a specific time t, the following steps are undertaken: (a) The interface position X 1 (t) (or X 2 (t)) is computed using Equation 14 (or Equation 16).(b) The velocity V 1 (t) (or V 2 (t)) is calculated based on Equation 15 (or Equation 17).(c) Assuming P out = 0, the input pressure P ent is determined using Equation A2.(d) The interface pressure P * 1 (or P * 2 ) is derived from Equation A1.Thus, at time t, the pressure within fracture i diminishes linearly in the initial segment (0 ≤ x ≤ X i (t)), ranging from P ent at x = 0 to P * i at x = X i (t).Subsequently, in the second segment (X i (t) ≤ x ≤ L), it continues to decrease linearly from P * i at x = X i (t) to P out = 0 at (x = L).
In this initial phase, since neither contaminant front has reached the exit point, the outlet concentration is consequently as follows:

The Flow Regime During the Second Phase (t 1 < t ≤ t 2 )
During the second phase, the contaminant from F1 has arrived at the outlet, while in F2, it has yet to reach the outlet, a situation assumed to occur at t = t 2 .Under these circumstances, Equation 13, applicable to both fractures, is expressed as follows: Water Resources Research As for the first phase, detailed calculations given in Appendix B (see Equation B8) yield: The velocity V 2 (t) is calculated as follows: The velocity in F1 varies with time and is obtained by employing Equation 5 as follows: The time t 2 when the contaminant in F2 reaches the outlet (X 2 (t 2 ) = L) is estimated by: In this second phase, the pressure distribution at a given time t in F1 linearly decreases from P ent at (x = 0) to P out = 0 at (x = L).The pressure P ent changes over time, and its value at a specific time t can be determined by applying Equations 26 and 23.In the first segment (0 ≤ x ≤ X i (t)) of F2, which is filled with the injected fluid, the pressure drops linearly from P ent at (x = 0) to P * 2 at (x = X 2 (t)).The interface position X 2 (t) is calculated using Equation 24, and the pressure at this interface, P * 2 , is derived from Equation 13.Following this, in the latter segment (X i (t) ≤ x ≤ L) of F2, the pressure at time t diminishes from P * 2 at (x = X 2 (t)) to P out = 0 at (x = L).
During the second phase, the fluid leaving F1 has a concentration of 1, whereas the fluid leaving F2 has a zero concentration, the mixing of both fluids creates an outlet concentration that varies with time:

The Flow Regime During the Third Phase (t ≥ t 2 )
In this third phase, both fractures are filled by the contaminant.In this case, Equation 13 becomes: As for the tracer case, we obtain fixed velocities using Equation 7.
In this final phase, the pressure remains constant over time and exhibits a linear variation from P ent at x = 0 to P out at x = L. Assuming P out = 0, the input pressure P ent is higher than for the tracer test case attributable to the reduction in mobility resulting from the increased viscosity.P ent is calculated using Equation 7 as: Water Resources Research 10.1029/2023WR035741 Furthermore, the fluid arriving from both fractures has a concentration of 1, as a consequence, the outlet concentration is:

Numerical Experiments
The test problem was addressed using two different methodologies: an analytical approach with the solution we formulated, and numerical techniques involving two distinct 2D in-house codes.The first numerical approach utilized the advanced Discontinuous Galerkin Finite Element (DGFE) method, as developed by Younes et al. (2022), which integrates highorder time integration strategies through the Method of Lines (MOL).The second approach employed the Finite Volume (FV) method.Implementing these methods enabled a thorough comparative analysis of different solution techniques.In the numerical simulations, the discretization of both the matrix and fracture continua was achieved using triangular elements on an unstructured mesh.This mesh was refined around the fractures to account for their narrow aperture, comprising a total of 12,998 triangles, as shown in Figure 1.The fracture aperture ratio is set at δ = 0.5, reflecting the relative size of the aperture, and the viscosity ratio is fixed at β = 0.138, indicating the comparative viscosity between the injected and displaced fluids.Additional input parameters are listed in Table 1.

Results and Discussion
In the subsequent subsections, we present the analytical solutions formulated for the specific problem being examined.These solutions are then compared with numerical solutions, highlighting the agreements and differences, as well as discussing the potential strengths and weaknesses of each method.Moreover, we concentrate on analyzing the BTC data to discern and estimate the pertinent parameters of the system under investigation.

Analytical Solutions
For the tracer test case, the analytical solution (Equations 7-12) yield fixed values V 1 = 4.44 × 10 4 m/s, V 2 = 1.11 × 10 4 m/s, t 1 = 1,170 s, t 2 = 4,680 s and C L (t)| t 1 ≤ t ≤ t 2 = 0.889.For the viscous test case, Equation 18yields t 1 = 1,266.93s indicating that the fluid in F1 experiences a delay, arriving approximately 100 s later at the outlet compared to the tracer test case.To determine the arrival time of the injected fluid in F2 at the outlet, Equation 27 is solved, resulting in t 2 = 3,334.88s which signifies that the fluid in F2 moves significantly faster, arriving around 1,300 s earlier at the outlet compared to the tracer test case.The analytical expressions for the velocity in the first fracture V 1 (t) are given by Equation 15for the first phase (t ≤ t 1 ), Equation 26for the second phase (t 1 ≤ t ≤ t 2 ), and Equation 29for the third phase (t 2 ≤ t).Similarly, the analytical expressions for the velocity in the second fracture are provided by Equation 17, Equation 25, and Equation 29for the three successive phases, respectively.
The pressure evolution at the middle of F1 and F2 are plotted in Figure 3 for the tracer and the viscous test cases.For the tracer test case, the same mid-fracture pressure is observed in F1 and F2.The mid-fracture pressure is constant during time and corresponds to 1 2 P ent where the pressure at the input, calculated using Equation 8, yields P ent = 0.69 [Pa].For the viscous test case, the mid-fracture F1 has a very small pressure decrease until it is reached by the injected fluid (Figure 3).Subsequently, the pressure rises markedly, following an almost linear trajectory, up until t 1 , at which point F1 becomes completely occupied by the injected fluid.Following this, the pressure at the midpoint of fracture F1 shows only a slight change until t 2 , the point at which both F1 and F2 are fully occupied by the injected fluid.Conversely, at the outset, the pressure at the midpoint of fracture F2 climbs moderately until the injected fluid reaches this midpoint.The pressure then experiences a notable rise until t 2 .Ultimately, once both fractures are  entirely filled with the injected fluid, the pressures at their midpoints converge, stabilizing at a constant value equal to 1 2 P ent , where P ent is determined using Equation 30 and is found to be 5.03 [Pa].Therefore, the final entry pressure in the case of higher viscosity is 7.26 times greater than that in the tracer case, reflecting the viscosity ratio between the injected and displaced fluids.

Comparison With Numerical Solutions
Figures 4a and 4b display the changes in velocity within fractures F1 and F2, respectively, as simulated through both numerical and analytical methods.At the experiment's onset, the fluid velocity in the wider fracture F1 is four times faster compared to the narrower fracture F2.This difference is attributed to the fact that F1's aperture is double the size of F2's.Consequently, the permeability (and similarly, the mobility) of F1 is quadruple that of F2, corresponding to the square of the aperture ratio.
When the viscous fluid progresses, the velocity in F1 gradually decreases, while the velocity in F 2 concurrently increases.In fact, as X 1 (t) increases, the equivalent mobility in F1, K eq 1 = , decreases which induces a reduction of the velocity inside F1.Consequently, the arrival time of the injected fluid in F1 at the outlet is t 1 = 1,266.93s which is slower than for the tracer test case.At the same time, because of the mass conservation of the total flow rate, an increase of the velocity in F2 is observed.Once, F1 is entirely filled by the injected fluid, it's equivalent mobility becomes constant.When the injected fluid progresses in F2, the equivalent permeability K eq 2 decreases which induces a reduction of the velocity inside F2 and by compensation, the velocity in F1 begins to rise, until the injected fluid in F2 reaches the outlet at the time t 2 = 3,334.88s.Subsequently, both fractures are filled with the same fluid (the injected fluid), resulting in the velocity in F 1 once again being four times higher than that in F 2 .A remarkable agreement is observed between the DGFE and the analytical solutions (Figure 4).However, a slight discrepancy arises when comparing the solution derived from a numerical model employing the FV method with the analytical solution.
Note that during the first phase (t ≤ t 1 ), the amount of reduction of the flow rate inside F1, causing a decrease of V 1 , produces, by compensation, the same amount of flow rate rise in F2, which induces a more significant enhance of V 2 because of the difference of the aperture of the fractures.Consequently, when compared to the tracer test case, the advance of t 2 is much more significant than the delay of t 1 .This is an important feature of injecting viscous fluids into fractured media, as they accelerate the crossing of small-aperture fractures.This feature is beneficial, for example, when injecting viscous surfactants to improve hydrocarbon recovery in a fractured reservoir or for decontamination of fractured soils.
Figure 5 displays the analytical and numerical breakthrough curves.For the viscous case represented in Figure 5b, the outlet concentration initially remains at zero until the injected fluid from F1 reaches the outlet at the time t 1 = 1,266.93s.At this point, a distinct and sharp front is observed.Similarly, a second sharp front is observed when the injected fluid from F2 reaches the outlet at the time t 2 = 3,334.88s.Following these breakthrough events, the normalized outlet concentration reaches a value of one.However, during the intermediate period between t 1 and t 2 , the outlet concentration experiences variations.This is due to the mixing of the fluid from F1, which has a concentration of 1 and a flux of e 1 V 1 (t) , with the fluid from F2, which has a concentration of 0 and a flux of e 2 V 2 (t) .As a result, the outlet concentration is not constant between t 1 and t 2 because the velocity in the fractures changes over time.This variation in concentration explains why the plateau observed in the tracer test case, as depicted in Figure 5a between the two sharp fronts is not observed in the viscous case represented in Figure 5b.Because of the negligence of the diffusion process, the analytical solution in Figure 5 exhibits sharp arriving fronts for both tracer and viscous cases.These fronts are very well reproduced by the DGFE numerical model, while the FV model fails to adequately reproduce these fronts and adds artificial numerical diffusion smearing the sharp BTCs.

Parameter Identifiability From the Breakthrough Curve
BTCs reflect flow behavior and exhibit sensitivity to various model parameters, such as fracture aperture, fluid viscosity, and transport properties.Consequently, BTCs have demonstrated their utility in estimating flow parameters within porous media (Hoffmann et al., 2022;Tran and Jha, 2020) and characterizing fracture properties (Hyman & Dentz, 2021;Srinivasan et al., 2019).In this sub-section, our focus lies in assessing the identifiability of the following parameters from BTCs: the aperture of F1 (e 1 ), the aperture of F2 (e 2 ), and the viscosity ratio between the injected and displaced fluids (η = μ 1 /μ 0 ).We aim to accomplish this using noisy measurements of the BTCs while assuming equal fluid densities for the injected and displaced fluids.
The vector of unknown parameters is represented by ξ, and we perform parameter estimation using a Bayesian approach by combining prior parameter information with the observed BTCs.This process simulates many samples with the developed analytical model and allows us to determine the posterior probability distribution functions (PDFs) of the model parameters.
To accomplish this, we rely on the widely utilized Markov chain Monte Carlo technique (MCMC), which has been successfully employed by several authors in the field of Hydrogeology (e.g., Du et al., 2022;Linde et al., 2017;Rajabi & Ataie-Ashtiani, 2016;Wang et al., 2019;Wei et al., 2023;Younes et al., 2016).MCMC is a Bayesian inference method where parameters are treated as random variables.MCMC generates random sequences of parameter sets that gradually converge toward the desired target distribution.By obtaining statistical measures such as the mean and standard deviation from these distributions, we can estimate the mean parameter values and determine their confidence intervals.This allows us to effectively characterize parameter uncertainty.Applying the Bayes theorem, we express the posterior density function of the parameters conditioned on the observations as follows (Turkman et al., 2019): In which p (ξ|y mes ) is the likelihood function measuring how well the model outputs are in agreement with the observations y mes , and p(ξ) is the prior PDF of ξ, which encapsulates any prior knowledge about the unknown parameters.
In this study, we assume that the prior distributions for the three unknown parameters are independent of each other and follow uniform distributions.To account for a wide range of potential parameter values, large prior intervals are chosen for all parameters.We also assume that the measurement errors are normally distributed and independent from each other.To investigate the identifiability of the parameters, we consider two scenarios with different levels of measurement errors.The first scenario involves high measurement errors, where the standard deviation (σ e ) is equal to 0.05.The second scenario corresponds to moderate measurement errors, where σ e = 0.02.These error levels are chosen specifically for normalized concentrations bounded between 0 and 1, which are the range of interest in this study.
Setting the calibration problem in a Bayesian framework yields the following posterior PDF (Turkman et al., 2019): where 2 is the sum of the squares of differences between the observed (C (k) mes ) and predicted (C (k)  mod ) concentrations at the time t k , and N c is the overall number of observed concentrations.The MCMC sampler generates a new candidate ξ i from a proposal distribution q (ξ i |ξ i 1 ) which only depends on the previously accepted candidate.The new candidate can be accepted or rejected depending on the Hasting ratio α H (Yildirim, 2012): In this work, we use the DREAM (ZS) MCMC sampler, based on the Metropolis-Hastings algorithm (Vrugt et al., 2003).This sampler runs multiple different Markov chains in parallel and uses a discrete proposal distribution to evolve the sampler to the posterior distribution.DREAM (ZS) has excellent performance on complex, multimodal search problems (Laloy & Vrugt, 2012).The sampler is run with three parallel chains and the results are considered stationary if the chains are not auto-correlated and if the Gelman and Rubin (1992) criterion is verified (R stat ≤1.2).The set of parameters corresponding to the maximum a posteriori (MAP) value is then defined as (Turkman et al., 2019): The perturbed observations (y mes ) used for model calibration as conditioning information are illustrated in Figure 6.These observations represent highly noised measurements, with a standard deviation equal to 0.05.The vector consists of 120 measurements of the outlet concentration recorded over time.
To accommodate both error measurements σ e = 0.02, 0.05, the MCMC sampler is configured with three parallel chains, each comprising a total of 6,000 runs.In the subsequent analysis, the last 25% of the runs, which adequately capture the model-observation fit, are utilized to estimate the joint posterior distribution.
The MCMC sampler achieves convergence after approximately 3,000 model runs for both σ e values.The estimated mean values and the corresponding 95% Confidence Intervals (CIs) for each parameter are presented in Table 2.Note that the CI is computed based on the standard deviation, assuming a Gaussian posterior distribution.
The results demonstrate that the apertures of the two fractures are well estimated for both measurement errors (σ e = 0.05 and σ e = 0.02).The 95% CIs associated with the aperture values are small, indicating high confidence in the estimates.It is noteworthy that as the measurement errors decrease, the uncertainty in the estimated aperture values also decreases, aligning with expectations.On the other hand, the estimation of the viscosity ratio is comparatively less accurate than that of the fracture apertures.When confronted with larger measurement errors (σ e = 0.05), the uncertainty associated with the viscosity ratio (η) is 40% larger than that observed with smaller measurement errors (σ e = 0.02).
The MCMC results for the scenario with small measurement errors are presented in Figure 7.The diagonal plots display the inferred posterior parameter distributions, revealing relative bell-shaped distributions for all parameters.The off-diagonal scatterplots illustrate the pairwise correlations observed in the MCMC samples.It is worth noting that there are negligible correlations between most of the parameters.However, a moderate correlation (r = 0.94) is observed between e 2 and η, indicating a relationship between these two parameters.

Conclusion
In this study, we formulated a novel analytical solution that accounts for VVF in porous media with explicitly represented fractures.The domain of interest is a rock matrix with notably low permeability and considered as inactive, bisected by two fractures with distinct apertures.We analyzed the scenario where the viscosities of the injected and displaced fluids vary, leading to temporal changes in fluid velocities.This results in three unique phases during the infiltration process, for which we develop the analytical solution separately.A comparison of our analytical solution with the DGFE numerical solution indicates a high level of consistency in terms of general flow behavior, velocity, and BTC.However, discrepancies emerge when comparing our solution with the FV  method.Despite this, our analytical model effectively approximates VVF in fractured rocks while drastically reducing the computational time relative to numerical models.Moreover, we evaluate the identifiability of parameters such as fracture apertures and viscosity ratios between injected and displaced fluids.This evaluation, informed by BTC from our analytical solution, employs the MCMC technique with the DREAM (ZS) sampler.The apertures of the two fractures are estimated with high confidence for varying measurement errors, with reduced error correlating with a lower uncertainty in the estimated values.Conversely, viscosity ratio estimation is less precise than aperture estimation, with larger measurement errors leading to greater uncertainty.This finding underscores the significance of accurate measurements for minimizing parameter uncertainties.

Appendix B
The velocity equations given in Equation 23 yield: Combining Equation B1 with Equation 5gives: Using the change of variables: We obtain: The integration of this equation leads to: which yields: where cte 1 is obtained using the initial condition X 2 (t 1 ) = X 2 * which yields: and finally, we obtain: )

Figure 1 .
Figure 1.Conceptual framework of the test case and 2D computational mesh.

Figure 2 .
Figure 2. Schematic of the flow regime during the first phase.

Figure 3 .
Figure 3. Mid-fracture pressure evolution in F1 and F2 for the tracer and the viscous test cases.

Figure
Figure Comparison between the velocities obtained through analytical calculations and those estimated numerically in F1 (a) and F2 (b).

Figure 5 .
Figure 5.Comparison of the BTCs obtained through analytical calculations and those estimated numerically for the tracer (a) and viscous (b) test cases.

Figure 6 .
Figure 6.Noised BTC with a standard deviation of 0.05.

Table 1
Flow and Transport Parameters for the Test Case Problem

Table 2
Prior Intervals, Estimated Mean Values, and Confidence Intervals (CIs)for the Three Unknown Parameters YOUNES ET AL.