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

Anis Younes, Mohammad Mahdi Rajabi, Fatemeh Rezaiezadeh Roukerd, Marwan Fahs

► To cite this version:

Anis Younes, Mohammad Mahdi Rajabi, Fatemeh Rezaiezadeh Roukerd, Marwan Fahs. An Analytical Solution for Variable Viscosity Flow in Fractured Media: Development and Comparative Analysis With Numerical Simulations. Water Resources Research, 2024, 60, 10.1029/2023WR035741. insu-04725755

HAL Id: insu-04725755

https://insu.hal.science/insu-04725755v1

Submitted on 8 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Water Resources Research®

## RESEARCH ARTICLE

10.1029/2023WR035741

**Key Points:**

* Analytical model developed for variable viscosity flow in fractured media, overcoming the common uniform viscosity assumption

* Analytical and numerical solutions demonstrate overall consistency, with minor discrepancies observed

* Breakthrough curves from the analytical solution, using Markov Chain Monte Carlo, enable confident estimation of fracture apertures

**Correspondence to:**

M. M. Rajabi,
mmrajabi@modares.ac.ir

**Citation:**

Younes, A., Rajabi, M. M., Rezaiezadeh Roukerd, F., & Fahs, M. (2024). An analytical solution for variable viscosity flow in fractured media: Development and comparative analysis with numerical simulations. *Water Resources Research*, 60, e2023WR035741. https://doi.org/10.1029/2023WR035741

Received 20 JULY 2023

Accepted 26 FEB 2024

**Author Contributions:**

**Conceptualization:** Anis Younes, Mohammad Mahdi Rajabi

**Formal analysis:** Anis Younes

**Methodology:** Anis Younes

**Writing – original draft:** Anis Younes, Mohammad Mahdi Rajabi

**Writing – review & editing:** Fatemeh Rezaiezadeh Roukerd, Marwan Fahs

© 2024. The Authors.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_2.png" /></figure>

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

Anis Younes¹ ID, Mohammad Mahdi Rajabi²,³ ID, Fatemeh Rezaiezadeh Roukerd², and Marwan Fahs¹ ID

¹ITES, University of Strasbourg, CNRS, ENGEES, Strasbourg, France, ²Faculty of Civil and Environmental Engineering, Tarbiat Modares University, Tehran, Iran, ³Department of Engineering, Institute of Computational Engineering, University of Luxembourg, Esch Sur Alzette, Luxembourg

**Abstract** 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.

### 1. 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 (Anima-sau, 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 (Ananhaswamy & 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).

VVF processes have been largely investigated in unfractured domains (Zidane, 2023; Graf & Therrien, 2007; Shafabakhsh et al., 2021). Past studies have examined the behavior of VVF in fractured domains with respect to CO₂ storage (Kim et al., 2019; Raad & Hassanzadeh, 2018; Shafabakhsh et al., 2021; Wang et al., 2022), saltwater intrusion (Grillo et al., 2010; Hosseini et al., 2020; Koohbor et al., 2020; Sebben et al., 2015; Werner et al., 2013), nuclear waste disposal (Follin & Stigsson, 2014), thermal convection (Gruais & Poliševski, 2021; Mezon et al., 2018), and oil and gas production (Shen et al., 2016).

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.

## 2. 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, multi-continuum 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), and Hyman 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
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 and Wang (2023), and Toller (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.

### 3. 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):

$$ \nabla \cdot q = 0, \qquad (1) $$

$$ q = -\frac{k}{\mu} \nabla P, \qquad (2) $$

where P is the pressure [Pa], q is the Darcy velocity [LT⁻¹], μ is the fluid dynamic viscosity [ML⁻¹T⁻¹], k is the permeability of the fracture [L²], 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:

$$ \frac{\partial C}{\partial t} + V \nabla C = 0, \qquad (3) $$

where V = q/θ is the fluid velocity [LT⁻¹], θ 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):

$$ k = \frac{e^2}{12}, \qquad (4) $$

which is often referred to as the cubic law (Witherspoon et al., 1980).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Conceptual framework of the test case and 2D computational mesh.</figcaption></figure>

# 4. 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.

## 4.1. 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_{\text{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 $\mu_1$, whereas the displaced fluid has a viscosity of $\mu_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.

## 4.2. 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^2T^{-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:

$$Q = Q_1 + Q_2 = e_1 V_1 + e_2 V_2, \qquad (5)$$

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 $\mu_0$, implies that the pressure has a linear variation along each fracture. Assuming the pressure difference between the entry and outlet $\Delta P = P_{\text{ent}} - P_{\text{out}}$ [Pa], Equation 2 gives:

$$\frac{\Delta P}{\mu_0 L} = \frac{Q_1}{e_1 k_1} = \frac{Q_2}{e_2 k_2}, \qquad (6)$$

Using Equation 5 and assuming that $\delta = e_2/e_1$, we obtain:

$$V_1 = \frac{Q}{e_1(1 + \delta^3)}, \qquad (7)$$

$$V_2 = \frac{Q}{e_2(1 + \delta^3)},$$

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_{\text{ent}}$ at $x = 0$ and transitioning to $P_{\text{out}}$ at $x = L$. Assuming $P_{\text{out}} = 0$, the input pressure $P_{\text{ent}}$ can be calculated using Equation 6 as follows:

$$P_{\text{ent}} = \frac{Q\mu_0 L}{e_1 k_1 (1 + \delta^3)}, \qquad (8)$$

The advection Equation 3 implies that the concentration $C$ is constant along the characteristic $x(t) = X(t; \tilde{x}, \tilde{t})$ passing through the point $(\tilde{x}, \tilde{t})$ and determined by the following system of equations:

$$\left\{ \begin{array}{l} \frac{dx(t)}{dt} = V(x(t),t), \\ x(\tilde{t}) = \tilde{x}, \end{array} \right. \qquad (9)$$

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_{\text{inj}}$ in a domain initially free of contaminant, the concentration in the fracture is either $C(x,t) = 0$ if $(t^*<0)$ or $C(x,t) = C_{\text{inj}}$ if $(t^*\ge0)$.

The mass conservation of the contaminant at the outlet reflects the ideal mixing of the fluids arriving from F1 and F2:

$$QC_L(t) = Q_1 C_1(L,t) + Q_2 C_2(L,t), \qquad (10)$$

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:

$$C_L(t) = \frac{C_1(L,t)}{(1 + \delta^3)} + \frac{C_2(L,t)}{(1 + \delta^3)}, \qquad (11)$$

The breakthrough curve (BTC) illustrates the arrival of a rapid front from fracture 1 (F1) at the outlet at $t_1 = e_1(1 + \delta^3) (L/Q)$ followed by a slower front from fracture 2 (F2) at $t_2 = e_2(1 + \delta^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,
<figure><img src="image_4.png" /><figcaption>Figure 2. Schematic of the flow regime during the first phase.</figcaption></figure>

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:

$$C_L(t)_{t_1 < t < t_2} = \frac{1}{(1 + \delta^3)}, \qquad (12)$$

### 4.3. 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.

#### 4.3.1. 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 $\mu_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.

#### 4.3.2. The Flow Regime During the First Phase ($0 < t \le t_1$)

In the first phase, fracture i, which is L in length, contains the injected fluid with a viscosity of $\mu_1$ in its segment $X_i(t)$ and the displaced fluid, having a viscosity of $\mu_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 \le x \le X_i(t)$) of the fracture, which is filled with the injected fluid, and similarly linear in the latter segment ($X_i(t) \le x \le L$), occupied by the displaced fluid. The continuity of flux across these two segments within fracture i is expressed as follows:

$$Q_i(t) = e_i K_i^1 \frac{P_{\text{ent}} - P_i^*}{X_i(t)} = e_i K_i^2 \frac{P_i^* - P_{\text{out}}}{(L - X_i(t))}, \qquad (13)$$

with the motilities $K_i^1 = \frac{k_i}{\mu_1} = \frac{(e_i)^2}{12\mu_1}$ and $K_i^2 = \frac{k_i}{\mu_0} = \frac{(e_i)^2}{12\mu_0}$.

Detailed calculations presented in Appendix A yield the position of the interface in F1, as detailed in Equation A14:

$$X_1(t) = -\sqrt{\lambda_2^1 t^2 + 2\lambda_1^1 t + a_0^2} - a_1 t + a_0, \qquad (14)$$

Therefore, the fluid velocity in F1 is:

$$V_1(t) = \frac{dX_1(t)}{dt} = \frac{\lambda_2^1 t + \lambda_1^1}{\sqrt{\lambda_2^1 t^2 + 2\lambda_1^1 t + a_0^2}} - a_1, \qquad (15)$$

The interface position in F2 is given by combining Equation 14 and Equation A5:
$$X_2(t) = \frac{Q}{e_2}t + \delta^{-1}\left(\sqrt{\lambda_2^1 t^2 + 2\lambda_1^1 t + a_0^2} + a_1 t - a_0\right), \quad (16)$$

The fluid velocity in F2 is obtained from Equation 5:

$$V_2(t) = \frac{Q}{e_2} \frac{V_1(t)}{\delta} = \frac{Q}{e_2} + \frac{1}{\delta} \left( \frac{\lambda_2^1 t + \lambda_1^1}{\sqrt{\lambda_2^1 t^2 + 2\lambda_1^1 t + a_0^2}} + a_1 \right), \quad (17)$$

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:

$$\sqrt{\lambda_2^1 t_1^2 + 2\lambda_1^1 t_1 + a_0^2} + a_1 t_1 - a_0 + L = 0, \quad (18)$$

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 14 and Equation 16 and 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:

$$X_2^* = \frac{Q}{e_2} \frac{L}{t_1 - \delta}, \quad (19)$$

At $t = t_1$, the velocity in the two fractures is:

$$V_1(t_1) = \frac{(\delta^4 - 1) a_1^2 t_1 + a_1 (b_0 - L)}{a_1 t_1 - a_0 + L}, \quad (20)$$

$$V_2(t_1) = \frac{Q (\delta^4 - 1) a_1^2 t_1 + a_1 (b_0 - L)}{e_2 (a_1 t_1 - a_0 + L)}, \quad (21)$$

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_{\text{ent}} = 0$, the input pressure $P_{\text{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 \le x \le X_i(t)$), ranging from $P_{\text{ent}}$ at $x = 0$ to $P_i^*$ at $x = X_i(t)$. Subsequently, in the second segment ($X_i(t) \le x \le L$), it continues to decrease linearly from $P_i^*$ at $x = X_i(t)$ to $P_{\text{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:

$$C_{\text{out}}(t)|_{t \le t_1} = 0, \quad (22)$$

### 4.3.3. The Flow Regime During the Second Phase ($t_1 < t \le 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:

$$V_1(t) = \frac{K_1^1}{L} (P_{\text{ent}} - P_{\text{out}}), \quad (23)$$
$$V_2(t) = \frac{K_2^1 K_2^2}{(K_2^1 (L - X_2(t))) + (K_2^2 X_2(t))} (P_{\text{ent}} - P_{\text{out}}),$$

As for the first phase, detailed calculations given in Appendix B (see Equation B8) yield:

$$X_2(t) = \frac{1}{(1-\beta)} \left( \sqrt{2\delta^2 \frac{Q}{e_1} L(1-\beta)(t-t_1) + ((1-\beta)X_2^* + \beta L + \delta^3 L)^2 - \beta L - \delta^3 L} \right), \quad (24)$$

The velocity $V_2(t)$ is calculated as follows:

$$V_2(t) = \frac{dX_2(t)}{dt} = \frac{\delta^2 \frac{Q}{e_1} L}{\sqrt{2\delta^2 \frac{Q}{e_1} L(1-\beta)(t-t_1) + ((1-\beta)X_2^* + \beta L + \delta^3 L)^2}}, \quad (25)$$

The velocity in F1 varies with time and is obtained by employing Equation 5 as follows:

$$V_1(t) = \frac{Q}{e_1} \left( 1 - \frac{\delta^3 L}{\sqrt{2\delta^2 \frac{Q}{e_1} L(1-\beta)(t-t_1) + ((1-\beta)X_2^* + \beta L + \delta^3 L)^2}} \right), \quad (26)$$

The time $t_2$ when the contaminant in F2 reaches the outlet ($X_2(t_2) = L$) is estimated by:

$$t_2 = t_1 + \frac{(1 + \delta^3)^2 L^2 - ((1-\beta)X_2^* + \beta L + \delta^3 L)^2}{2\delta^2 \frac{Q}{e_1} L(1-\beta)}, \quad (27)$$

In this second phase, the pressure distribution at a given time $t$ in F1 linearly decreases from $P_{\text{ent}}$ at ($x = 0$) to $P_{\text{out}} = 0$ at ($x = L$). The pressure $P_{\text{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 \le x \le X_1(t)$) of F2, which is filled with the injected fluid, the pressure drops linearly from $P_{\text{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_1(t) \le x \le L$) of F2, the pressure at time $t$ diminishes from $P_2^*$ at ($x = X_2(t)$) to $P_{\text{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:

$$C_{\text{out}}(t)|_{t_1 \le t \le t_2} = \frac{e_1 V_1(t)}{Q} = 1 - \frac{\delta^3 L}{\sqrt{2\delta^2 \frac{Q}{e_1} L(1-\beta)(t-t_1) + ((1-\beta)X_2^* + \beta L + \delta^3 L)^2}}, \quad (28)$$

#### 4.3.4. The Flow Regime During the Third Phase ($t \ge t_2$)

In this third phase, both fractures are filled by the contaminant. In this case, Equation 13 becomes:

$$\frac{V_2(t)}{V_1(t)} = \frac{K_2^1}{K_1^1} = \delta^2, \quad (29)$$

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_{\text{ent}}$ at $x = 0$ to $P_{\text{out}}$ at $x = L$. Assuming $P_{\text{out}} = 0$, the input pressure $P_{\text{ent}}$ is higher than for the tracer test case attributable to the reduction in mobility resulting from the increased viscosity. $P_{\text{ent}}$ is calculated using Equation 7 as:
<table><caption>**Table 1**
Flow and Transport Parameters for the Test Case Problem</caption><thead><tr><th>Parameters</th><th>Value</th></tr></thead><tbody><tr><td>Rock matrix permeability</td><td>k<sub>m</sub> = 10<sup>-15</sup> m<sup>2</sup></td></tr><tr><td>Storage coefficient</td><td>S = 10<sup>-10</sup>m<sup>-1</sup></td></tr><tr><td>Aperture of the left fracture</td><td>e<sub>1</sub> = 2 × 10<sup>-3</sup>m</td></tr><tr><td>Aperture of the right fracture</td><td>e<sub>2</sub> = 10<sup>-3</sup>m</td></tr><tr><td>Length of the fractures</td><td>L = 0.52 m</td></tr><tr><td>Injected flow rate</td><td>Q = 10<sup>-6</sup> m<sup>2</sup>/s</td></tr><tr><td>Dispersivities</td><td>α<sub>L</sub> = α<sub>r</sub> = 0m</td></tr><tr><td>Molecular diffusion</td><td>D<sub>m</sub> = 0 m<sup>2</sup>/s</td></tr><tr><td>Density of the displaced fluid</td><td>ρ<sub>0</sub> = 1000 Kg/m<sup>3</sup></td></tr><tr><td>Viscosity of the displaced fluid</td><td>μ<sub>0</sub> = 0.001 Kg/m.s</td></tr><tr><td>Density of the injected fluid</td><td>ρ<sub>1</sub> = 1000 Kg/m<sup>3</sup></td></tr><tr><td>Density of the injected fluid</td><td>μ<sub>1</sub> = 0.00726 Kg/m.s</td></tr></tbody></table>

$$P_{\text{ent}} = \frac{Q\mu_1 L}{e_1 k_1 (1 + \delta^3)}, \qquad (30)$$

Furthermore, the fluid arriving from both fractures has a concentration of 1, as a consequence, the outlet concentration is:

$$C_{\text{out}}(t)|_{t \ge t_2} = 1, \qquad (31)$$

## 5. 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 high-order 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 $\delta = 0.5$, reflecting the relative size of the aperture, and the viscosity ratio is fixed at $\beta = 0.138$, indicating the comparative viscosity between the injected and displaced fluids. Additional input parameters are listed in Table 1.

## 6. 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.

### 6.1. Analytical Solutions

For the tracer test case, the analytical solution (Equations 7–12) yield fixed values $V_1 = 4.44 \times 10^{-4}$ m/s, $V_2 = 1.11 \times 10^{-4}$ m/s, $t_1 = 1,170$ s, $t_2 = 4,680$ s and $C_L(t)|_{t_1 \le t \le t_2} = 0.889$. For the viscous test case, Equation 18 yields $t_1 = 1,266.93$ s 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.88$ s 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 15 for the first phase ($t \le t_1$), Equation 26 for the second phase ($t_1 \le t \le t_2$), and Equation 29 for the third phase ($t_2 \le t$). Similarly, the analytical expressions for the velocity in the second fracture are provided by Equation 17, Equation 25, and Equation 29 for 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 $\frac{1}{2}P_{\text{ent}}$ where the pressure at the input, calculated using Equation 8, yields $P_{\text{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
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Mid-fracture pressure evolution in F1 and F2 for the tracer and the viscous test cases.</figcaption></figure>

entirely filled with the injected fluid, the pressures at their midpoints converge, stabilizing at a constant value equal to $\frac{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.

## 6.2. 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 F2 concurrently increases. In fact, as $X_1(t)$ increases, the equivalent mobility in F1, $K_1^{eq} = \frac{e_1^2 L}{\mu_1 X_1(t) + \mu_0 (L - X_1(t))}$, 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.93$ s 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_2^{eq}$ 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.88$ s. Subsequently, both fractures are filled with the same fluid (the injected fluid), resulting in the velocity in F1 once again being four times higher than that in F2. 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 \le 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.93$ s. 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.88$ s. 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.
<figure><img src="image_6.png" /><figcaption>Figure 4. Comparison between the velocities obtained through analytical calculations and those estimated numerically in F1 (a) and F2 (b).</figcaption></figure>

## 6.3. 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 ($\eta = \mu_1/\mu_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 $\xi$, 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):

$$p(\xi | y_{\text{mes}}) \propto p(y_{\text{mes}} | \xi) p(\xi), \quad (32)$$

In which $p(\xi|y_{\text{mes}})$ is the likelihood function measuring how well the model outputs are in agreement with the observations $y_{\text{mes}}$, and $p(\xi)$ is the prior PDF of $\xi$, 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 ($\sigma_e$) is equal to 0.05. The second scenario corresponds to moderate measurement errors, where $\sigma_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):

$$p(\xi | y_{\text{mes}}) \propto \sigma_e^{-N_c} \exp \left( -\frac{SS_C(\xi)}{2\sigma_e^2} \right), \quad (33)$$

where $SS_C(\xi) = \sum_{k=1}^{N_c} (C_{\text{mes}}^{(k)} - C_{\text{mod}}^{(k)}(\xi))^2$ is the sum of the squares of differences between the observed ($C_{\text{mes}}^{(k)}$) and predicted ($C_{\text{mod}}^{(k)}$) concentrations at the time $t_k$, and $N_c$ is the overall number of observed concentrations. The
<figure><img src="image_7.png" /><figcaption>(b) Different viscosity of injected and displaced fluids.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison of the BTCs obtained through analytical calculations and those estimated numerically for the tracer (a) and viscous (b) test cases.</figcaption></figure>

MCMC sampler generates a new candidate $\xi^i$ from a proposal distribution $q(\xi^i|\xi^{i-1})$ which only depends on the previously accepted candidate. The new candidate can be accepted or rejected depending on the Hasting ratio $\alpha_H$ (Yildirim, 2012):

$$ \alpha_H = \min \left( 1, \frac{p(\xi^i | y_{\text{mes}}) q(\xi^i | \xi^{i-1})}{p(\xi^{i-1} | y_{\text{mes}}) q(\xi^{i-1} | \xi^i)} \right), \quad (34) $$

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
<figure><img src="image_9.png" /><figcaption>Figure 6. Noised BTC with a standard deviation of 0.05.</figcaption></figure>

are considered stationary if the chains are not auto-correlated and if the Gelman and Rubin (1992) criterion is verified ($R_{stat} \le 1.2$). The set of parameters corresponding to the maximum a posteriori (MAP) value is then defined as (Turkman et al., 2019):

$$ \xi^{\text{MAP}} = \arg \max(p(\xi|\mathbf{y}_{\text{mes}})), \quad (35) $$

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 $\sigma_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 $\sigma_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 ($\sigma_e = 0.05$ and $\sigma_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 ($\sigma_e = 0.05$), the uncertainty associated with the viscosity ratio ($\eta$) is 40% larger than that observed with smaller measurement errors ($\sigma_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 $\eta$, indicating a relationship between these two parameters.

## 7. 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

<table><caption>**Table 2**
Prior Intervals, Estimated Mean Values, and Confidence Intervals (CIs) for the Three Unknown Parameters</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="2">Prior intervals</th><th colspan="3">σ<sub>e</sub> = 0.05</th><th colspan="3">σ<sub>e</sub> = 0.02</th></tr><tr><th>Lower bound</th><th>Upper bound</th><th>Mean</th><th>Standard deviation</th><th>95% CI</th><th>Mean</th><th>Standard deviation</th><th>95% CI</th></tr></thead><tbody><tr><td>e<sub>1</sub> [cm]</td><td>0.05</td><td>1.0</td><td>0.2</td><td>1.4 10<sup>-3</sup></td><td>(0.197–0.203)</td><td>0.2</td><td>4.2 10<sup>-4</sup></td><td>(0.199–0.201)</td></tr><tr><td>e<sub>2</sub> [cm]</td><td>0.05</td><td>1.0</td><td>0.1</td><td>1.6 10<sup>-3</sup></td><td>(0.097–0.103)</td><td>0.1</td><td>8.9 10<sup>-4</sup></td><td>(0.098–0.102)</td></tr><tr><td>η [-]</td><td>1</td><td>10</td><td>7.27</td><td>1.68</td><td>(4–10)</td><td>7.57</td><td>1.2</td><td>(5.2–9.9)</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>Figure 7. MCMC solutions using BTC with measurements errors $\sigma_e = 0.02$. The diagonal plots represent the inferred posterior parameter distributions. The off-diagonal represents the pairwise correlations between parameters.</figcaption></figure>

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 A

From Equation 14, the pressure at the interface is:

$$P_i^* = \frac{\frac{K_i^1}{X_i(t)} P_{\text{ent}} + \frac{K_i^2}{(L-X_i(t))} P_{\text{out}}}{\left(\frac{K_i^1}{X_i(t)} + \frac{K_i^2}{(L-X_i(t))}\right)}, \quad (A1)$$

Substituting Equation A1 into Equation 14 gives the velocity $V_i(t)$ as:

$$V_i(t) = \frac{K_i^1 K_i^2}{(K_i^1 (L-X_i(t)) + K_i^2 X_i(t))} (P_{\text{ent}} - P_{\text{out}}), \quad (A2)$$

Writing Equation A2 for $i = 1, 2$ and denoting $\frac{K_1^1}{K_1^2} = \frac{K_2^1}{K_2^2} = \frac{\mu_0}{\mu_1} = \beta$ gives:

$$\frac{V_2(t)}{V_1(t)} = \frac{\delta^2 \beta (L - X_1(t)) + X_1(t)}{\beta (L - X_2(t)) + X_2(t)}, \quad (A3)$$
Substituting Equation A3 into Equation 5 gives:

$$V_1(t) = \frac{dX_1(t)}{dt} = \frac{Q/e_1}{1 + \delta^3 \frac{(1-\beta)}{(1-\beta)} \frac{X_1(t)+\beta L}{X_2(t)+\beta L}}, \quad (A4)$$

Integrating Equation 5 and using the initial conditions $X_1(0) = X_2(0) = 0$ gives:

$$e_1 X_1(t) + e_2 X_2(t) = Q t, \quad (A5)$$

Substituting Equation A5 into Equation A4 gives:

$$\frac{dX_1(t)}{dt} = \frac{-\delta^{-1} \frac{Q(1-\beta)}{e_1} X_1(t) + \frac{Q^2}{e_1 e_2} (1-\beta)t + \beta \frac{Q}{e_1} L}{(\delta^3 - \delta^{-1})(1-\beta)X_1(t) + \frac{Q(1-\beta)}{e_2} t + \beta L(1+\delta^3)}, \quad (A6)$$

Which writes:

$$\frac{dX_1(t)}{dt} = -a_1 \frac{(X_1(t) + a_1 t - a_0) - a_1 \delta^4 t + a_0 - b_0}{(X_1(t) + a_1 t - a_0)}, \quad (A7)$$

where $a_0 = \frac{(1+\delta^3)}{(\delta^3 - \delta^{-1})} \frac{\beta L}{(\beta-1)}$, $a_1 = \frac{Q}{(\delta^4 - 1) e_1}$, and $b_0 = \frac{\delta \beta L}{(1-\beta)}$. Using the following change of variables:

$$Y_1(t) = X_1(t) + a_1 t - a_0, \quad (A8)$$

Equation A7 becomes:

$$\frac{dY_1(t)}{dt} - a_1 = -a_1 \frac{Y_1(t) - a_1 \delta^4 t + a_0 - b_0}{Y_1(t)}, \quad (A9)$$

which simplifies to:

$$Y_1(t) \frac{dY_1(t)}{dt} = \lambda_2^1 t + \lambda_1^1, \quad (A10)$$

where $\lambda_1^1 = a_1(b_0-a_0)$ and $\lambda_2^1 = \delta^4 a_1^2$. The integration of Equation A10 yields:

$$Y_1(t) = -\sqrt{\lambda_2^1 t^2 + 2\lambda_1^1 t + \lambda_0^1}, \quad (A11)$$

where $\lambda_0^1$ is a constant. Thus, Equation A8 gives:

$$X_1(t) = -\sqrt{\lambda_2^1 t^2 + 2\lambda_1^1 t + \lambda_0^1 - a_1 t + a_0}, \quad (A12)$$

Considering the initial condition $X_1(0) = 0$, we obtain:

$$\lambda_0^1 = a_0^2, \quad (A13)$$

and Equation A12 becomes:

$$X_1(t) = -\sqrt{\lambda_2^1 t^2 + 2\lambda_1^1 t + a_0^2 - a_1 t + a_0}, \quad (A14)$$
## Appendix B

The velocity equations given in Equation 23 yield:

$$ \frac{V_2(t)}{V_1(t)} = \frac{K_2^1 K_2^2 L}{K_1^1 (K_2^1 (L-X_2(t)) + K_2^2 X_2(t))}, \quad (B1) $$

Combining Equation B1 with Equation 5 gives:

$$ V_2(t) = \frac{Q \delta^2 \frac{L}{e_1}}{(1-\beta) X_2(t) + \beta L + \delta^3 L}, \quad (B2) $$

Using the change of variables:

$$ Y_2(t) = (1-\beta)X_2(t) + \beta L + \delta^3 L, \quad (B3) $$

We obtain:

$$ Y_2(t) \frac{dY_2(t)}{dt} = \delta^2 \frac{Q}{e_1} L(1-\beta), \quad (B4) $$

The integration of this equation leads to:

$$ Y_2(t) = \sqrt{2\delta^2 \frac{Q}{e_1} L(1-\beta)t + cte_1}, \quad (B5) $$

which yields:

$$ X_2(t) = \frac{1}{(1-\beta)} \left( \sqrt{2\delta^2 \frac{Q}{e_1} L(1-\beta)t + cte_1} - \beta L - \delta^3 L \right), \quad (B6) $$

where $cte_1$ is obtained using the initial condition $X_2(t_1) = X_2^*$ which yields:

$$ cte_1 = ((1-\beta)X_2^* + \beta L + \delta^3 L)^2 - 2\delta^2 \frac{Q}{e_1} L(1-\beta)t_1, \quad (B7) $$

and finally, we obtain:

$$ X_2(t) = \frac{1}{(1-\beta)} \left( \sqrt{2\delta^2 \frac{Q}{e_1} L(1-\beta)(t-t_1) + ((1-\beta)X_2^* + \beta L + \delta^3 L)^2 - \beta L - \delta^3 L} \right). \quad (B8) $$

## Data Availability Statement

Data were not used, nor created for this research.

## Acknowledgments

This work was supported by the French National Center of Research (CNRS) and Strasbourg University.

## References

Ackerer, P., Younes, A., & Mose, R. (1999). Modeling variable density flow and solute transport in porous medium: 1. Numerical model and verification. *Transport in Porous Media*, 35(3), 345–373. https://doi.org/10.1023/A:1006564309167

Adler, P. M., Thovert, J. F., & Mourzenko, V. V. (2013). *Fractured porous media*. Oxford University Press.

Agrawal, P., Dadheech, P. K., Jat, R. N., Bohra, M., Nisar, K. S., & Khan, I. (2020). Lie similarity analysis of MHD flow past a stretching surface embedded in porous medium along with imposed heat source/sink and variable viscosity. *Journal of Materials Research and Technology*, 9(5), 10045–10053. https://doi.org/10.1016/j.jmrt.2020.07.023

Ahmed, M. I., Abd-Elmegeed, M. A., & Hassan, A. E. (2019). Modelling transport in fractured media using the fracture continuum approach. *Arabian Journal of Geosciences*, 12(5), 172. https://doi.org/10.1007/s12517-019-4314-3
Al-Khafajy, D. G. (2016). Effects of heat transfer on MHD oscillatory flow of Jeffrey fluid with variable viscosity through porous medium. *Advances in Applied Science Research*, 7(3), 179-186.

Ananthaswamy, V., & Uma Maheswari, S. (2015). Analytical expression for the hydrodynamic fluid flow through a porous medium. *International Journal of Automation and Control Engineering*, 4(2), 67-76.

Angot, P. (2021). Solvability of the variable-viscosity fluid-porous flows coupled with an optimal stress jump interface condition.

Animasaun, I. L. (2015). Effects of thermophoresis, variable viscosity and thermal conductivity on free convective heat and mass transfer of non-darcian MHD dissipative Casson fluid flow with suction and nth order of chemical reaction. *Journal of the Nigerian Mathematical Society*, 34(1), 11-31.

Badea, L., Discacciati, M., & Quarteroni, A. (2010). Numerical analysis of the Navier-Stokes/Darcy coupling. *Numerische Mathematik*, 115(2), 195-227. https://doi.org/10.1007/s00211-009-0279-6

Baigereyev, D., Alimbekova, N., Berdyshev, A., & Madiyarov, M. (2021). Convergence analysis of a numerical method for a fractional model of fluid flow in fractured porous media. *Mathematics*, 9(18), 2179. https://doi.org/10.3390/math9182179

Berre, I., Doster, F., & Keilegavlen, E. (2019). Flow in fractured porous media: A review of conceptual models and discretization approaches. *Transport in Porous Media*, 130(1), 215-236. https://doi.org/10.1007/s11242-018-1171-6

Caucao, S., Gatica, G. N., Oyarzúa, R., & Šebestová, I. (2017). A fully-mixed finite element method for the Navier-Stokes/Darcy coupled problem with nonlinear viscosity. *Journal of Numerical Mathematics*, 25(2), 55-88. https://doi.org/10.1515/jnma-2015-0121

Chou, H. M., Wu, H. W., Lin, I. H., Yang, W. J., & Cheng, M. L. (2015). Effects of temperature-dependent viscosity on natural convection in porous media. *Numerical Heat Transfer, Part A: Applications*, 68(12), 1331-1350. https://doi.org/10.1080/10407782.2015.1012864

Choudhary, S., Devi, R., Mahajan, A., & Sunil (2023). Stability analysis in a couple-stress fluid layer with variable viscosity heated from below: Different conducting boundaries. *Chinese Journal of Physics*, 83, 94-102. https://doi.org/10.1016/j.cjph.2023.02.007

Cotta, R. M., Lisboa, K. M., & Zotin, J. L. Z. (2020). Integral transforms for flow and transport in discrete and continuum models of fractured heterogeneous porous media. *Advances in Water Resources*, 142, 103621. https://doi.org/10.1016/j.advwatres.2020.103621

Du, X., Du, C., Radolinski, J., Wang, Q., & Jian, J. (2022). Metropolis-hastings Markov chain Monte Carlo approach to simulate van genuchten model parameters for soil water retention curve. *Water*, 14(12), 1968. https://doi.org/10.3390/w14121968

Durlofsky, L. J. (1991). Numerical calculation of equivalent grid block permeability tensors for heterogeneous porous media. *Water Resources Research*, 27(5), 699-708. https://doi.org/10.1029/91wr00107

Ellahi, R., & Afzal, S. (2009). Effects of variable viscosity in a third grade fluid with porous medium: An analytic solution. *Communications in Nonlinear Science and Numerical Simulation*, 14(5), 2056-2072. https://doi.org/10.1016/j.cnsns.2008.05.006

Fahs, H., Hayek, M., Fahs, M., & Younes, A. (2014). An efficient numerical model for hydrodynamic parameterization in 2D fractured dual-porosity media. *Advances in Water Resources*, 63, 179-193. https://doi.org/10.1016/j.advwatres.2013.11.008

Fang, H., & Zhu, J. (2018). Simulation of groundwater exchange between an unconfined aquifer and a discrete fracture network with laminar and turbulent flows. *Journal of Hydrology*, 562, 468-476. https://doi.org/10.1016/j.jhydrol.2018.05.022

Follin, S., & Stigsson, M. (2014). A transmissivity model for deformation zones in fractured crystalline rock and its possible correlation to in situ stress at the proposed high-level nuclear waste repository site at Forsmark, Sweden. *Hydrogeology Journal*, 22(2), 299-311. https://doi.org/10.1007/s10040-013-1078-9

Gao, Q., & Ghassemi, A. (2020). Three-dimensional thermo-poroelastic modeling and analysis of flow, heat transport and deformation in fractured rock with applications to a lab-scale geothermal system. *Rock Mechanics and Rock Engineering*, 53(4), 1565-1586. https://doi.org/10.1007/s00603-019-01989-0

Gelman, A., & Rubin, D. B. (1992). [practical Markov chain Monte Carlo]: Rejoinder: Replication without conition. *Statistical Science*, 7(4), 503-511. https://doi.org/10.1214/ss/1177011148

Graf, T., & Boufadel, M. C. (2011). Effect of viscosity, capillarity and grid spacing on thermal variable-density flow. *Journal of Hydrology*, 400(1-2), 41-57. https://doi.org/10.1016/j.jhydrol.2011.01.025

Graf, T., & Simmons, C. T. (2009). Variable-density groundwater flow and solute transport in fractured rock: Applicability of the Tang et al. [1981] analytical solution. *Water Resources Research*, 45(2). https://doi.org/10.1029/2008WR007278

Graf, T., & Therrien, R. (2007). Coupled thermohaline groundwater flow and single-species reactive solute transport in fractured porous media. *Advances in Water Resources*, 30(4), 742-771. https://doi.org/10.1016/j.advwatres.2006.07.001

Grillo, A., Logashenko, D., Stichel, S., & Wittum, G. (2010). Simulation of density-driven flow in fractured porous media. *Advances in Water Resources*, 33(12), 1494-1507. https://doi.org/10.1016/j.advwatres.2010.08.004

Gruais, I., & Poliševski, D. (2021). Thermal flows in fractured porous media. *ESAIM: Mathematical Modelling and Numerical Analysis*, 55(3), 789-805. https://doi.org/10.1051/m2an/2020087

Hasona, W. M., El-Shekhipy, A. A., & Ibrahim, M. G. (2018). Combined effects of magnetohydrodynamic and temperature dependent viscosity on peristaltic flow of Jeffrey nanofluid through a porous medium: Applications to oil refinement. *International Journal of Heat and Mass Transfer*, 126, 700-714. https://doi.org/10.1016/j.ijheatmasstransfer.2018.05.087

Hayek, M. (2019). Accurate approximate semi-analytical solutions to the Boussinesq groundwater flow equation for recharging and discharging of horizontal unconfined aquifers. *Journal of Hydrology*, 570, 411-422. https://doi.org/10.1016/j.jhydrol.2018.12.057

Hirthe, E. M., & Graf, T. (2015). Fracture network optimization for simulating 2D variable-density flow and transport. *Advances in Water Resources*, 83, 364-375. https://doi.org/10.1016/j.advwatres.2015.07.001

Hoffmann, R., Maréchal, J. C., Selles, A., Dassargues, A., & Goderniaux, P. (2022). Heat tracing in a fractured aquifer with injection of hot and cold water. *Ground Water*, 60(2), 192-209. https://doi.org/10.1111/gwat.13138

Hosseini, N., Bajalan, Z., & Khoei, A. R. (2020). Numerical modeling of density-driven solute transport in fractured porous media with the extended finite element method. *Advances in Water Resources*, 136, 103453. https://doi.org/10.1016/j.advwatres.2019.103453

Huang, C. H., Zhu, H. Y., Wang, J. D., Han, J., Zhou, G. Q., & Tang, X. H. (2022). A FEM-DFN model for the interaction and propagation of multi-cluster fractures during variable fluid-viscosity injection in layered shale oil reservoir. *Petroleum Science*, 19(6), 2796-2809. https://doi.org/10.1016/j.petsci.2022.06.007

Hyman, J. D., & Dentz, M. (2021). Transport upscaling under flow heterogeneity and matrix-diffusion in three-dimensional discrete fracture networks. *Advances in Water Resources*, 155, 103994. https://doi.org/10.1016/j.advwatres.2021.103994

Hyman, J. D., Sweeney, M. R., Gable, C. W., Svyatsky, D., Lipnikov, K., & Moulton, J. D. (2022). Flow and transport in three-dimensional discrete fracture matrix models using mimetic finite difference on a conforming multi-dimensional mesh. *Journal of Computational Physics*, 466, 111396. https://doi.org/10.1016/j.jcp.2022.111396

Jayanthi, S., & Kumari, M. (2007). Effect of variable viscosity on non-Darcy free or mixed convection flow on a vertical surface in a non-Newtonian fluid saturated porous medium. *Applied Mathematics and Computation*, 186(2), 1643-1659. https://doi.org/10.1016/j.amc.2006.08.072
Jourde, H., Cornaton, F., Pistre, S., & Bidaux, P. (2002). Flow behavior in a dual fracture network. *Journal of Hydrology*, **266**(1-2), 99-119.
https://doi.org/10.1016/s0022-1694(02)00120-8

Kandasamy, R., Hashim, I., & Ruhaila (2008). Thermophoresis and chemical reaction effects on non-Darcy mixed convective heat and mass transfer past a porous wedge with variable viscosity in the presence of suction or injection. *Nuclear Engineering and Design*, **238**(10), 2699-2705. https://doi.org/10.1016/j.nucengdes.2008.05.010

Khoei, A. R., Saeedmonir, S., Hosseini, N., & Mousavi, S. M. (2023). An X-FEM technique for numerical simulation of variable-density flow in fractured porous media. *MethodsX*, **10**, 102137. https://doi.org/10.1016/j.mex.2023.102137

Kim, M., Kim, K. Y., Han, W. S., Oh, J., & Park, E. (2019). Density-driven convection in a fractured porous media: Implications for geological CO2 storage. *Water Resources Research*, **55**(7), 5852-5870. https://doi.org/10.1029/2019wr024822

Kim, M. C. K. (2014). Onset of buoyancy-driven convection in a variable viscosity liquid saturated in a porous medium. *Chemical Engineering Science*, **113**, 77-87. https://doi.org/10.1016/j.ces.2014.04.012

Koohbor, B., Fahs, M., Hoteit, H., Doummar, J., Younes, A., & Belfort, B. (2020). An advanced discrete fracture model for variably saturated flow in fractured porous media. *Advances in Water Resources*, **140**, 103602. https://doi.org/10.1016/j.advwatres.2020.103602

Kordilla, J., Sauter, M., Reimann, T., & Geyer, T. (2012). Simulation of saturated and unsaturated flow in karst systems at catchment scale using a double continuum approach. *Hydrology and Earth System Sciences*, **16**(10), 3909-3923. https://doi.org/10.5194/hess-16-3909-2012

Laloy, E., & Vrugt, J. A. (2012). High-dimensional posterior exploration of hydrologic models using multiple-try DREAM(ZS) and high-performance computing. *Water Resources Research*, **48**(1). https://doi.org/10.1029/2011WR010608

Linde, N., Ginsbourger, D., Irving, J., Nobile, F., & Doucet, A. (2017). On uncertainty quantification in hydrogeology and hydrogeophysics. *Advances in Water Resources*, **110**, 166-181. https://doi.org/10.1016/j.advwatres.2017.10.014

Liu, R., Li, B., Jiang, Y., & Huang, N. (2016). Mathematical expressions for estimating equivalent permeability of rock fracture networks. *Hydrogeology Journal*, **24**(7), 1623-1649. https://doi.org/10.1007/s10040-016-1441-8

Mezon, C., Mourzenko, V. V., Thovert, J. F., Antoine, R., Fontaine, F., Finizola, A., & Adler, P. M. (2018). Thermal convection in three-dimensional fractured porous media. *Physical Review A*, **97**(1), 013106. https://doi.org/10.1103/physreve.97.013106

O'Malley, D., Karra, S., Hyman, J. D., Viswanathan, H. S., & Srinivasan, G. (2018). Efficient Monte Carlo with graph-based subsurface flow and transport models. *Water Resources Research*, **54**(5), 3758-3766. https://doi.org/10.1029/2017wr022073

Prasad, V. R., Vasu, B., & Bég, O. A. (2011). Thermo-diffusion and diffusion-thermo effects on MHD free convection flow past a vertical porous plate embedded in a non-Darcian porous medium. *Chemical Engineering Journal*, **173**(2), 598-606. https://doi.org/10.1016/j.cej.2011.08.009

Raad, S. M. J., & Hassanzadeh, H. (2018). Onset of density-driven instabilities in fractured aquifers. *Physical Review A*, **97**(4), 043109. https://doi.org/10.1103/physreve.97.043109

Rajabi, M. M., & Ataie-Ashtiani, B. (2016). Efficient fuzzy Bayesian inference algorithms for incorporating expert knowledge in parameter estimation. *Journal of Hydrology*, **536**, 255-272. https://doi.org/10.1016/j.jhydrol.2016.02.029

Reddy, P. B. A., & Reddy, N. B. (2013). MHD Free Convection flow with variable viscosity and thermal diffusivity along a moving vertical plate embedded in a porous medium. *Acta Technica Corviniensis - Bulletin of Engineering*, **6**(3), 85.

Sabet, N., Hassanzadeh, H., & Abedi, J. (2018). A new insight into the stability of variable viscosity diffusive boundary layers in porous media under gravity field. *AIChE Journal*, **64**(3), 1083-1094. https://doi.org/10.1002/aic.15974

Schmid, K. S., Geiger, S., & Sorbie, K. S. (2013). Higher order FE-FV method on unstructured grids for transport and two-phase flow with variable viscosity in heterogeneous porous media. *Journal of Computational Physics*, **241**, 416-444. https://doi.org/10.1016/j.jcp.2012.12.017

Sebben, M. L., Werner, A. D., & Graf, T. (2015). Seawater intrusion in fractured coastal aquifers: A preliminary numerical investigation using a fractured henry problem. *Advances in Water Resources*, **85**, 93-108. https://doi.org/10.1016/j.advwatres.2015.09.013

Shafabakhsh, G., Sadeghnejad, M., & Ebrahimiia, R. (2021). Fracture resistance of asphalt mixtures under mixed-mode I/II loading at low-temperature: Without and with nano SiO2. *Construction and Building Materials*, **266**, 120954. https://doi.org/10.1016/j.conbuildmat.2020.120954

Shankar, B. M., Kumar, J., & Shivakumara, I. S. (2017). Stability of natural convection in a vertical layer of Brinkman porous medium. *Acta Mechanica*, **228**, 1-19. https://doi.org/10.1007/s00707-016-1690-6

Shen, W., Xu, Y., Li, X., Huang, W., & Gu, J. (2016). Numerical simulation of gas and water flow mechanism in hydraulically fractured shale gas reservoirs. *Journal of Natural Gas Science and Engineering*, **35**, 726-735. https://doi.org/10.1016/j.jngse.2016.08.078

Simacek, P., & Advani, S. G. (2003). A numerical model to predict fiber tow saturation during liquid composite molding. *Composites Science and Technology*, **63**(12), 1725-1736. https://doi.org/10.1016/s0266-3538(03)00155-6

Srinivasan, S., Karra, S., Hyman, J., Viswanathan, H., & Srinivasan, G. (2019). Model reduction for fractured porous media: A machine learning approach for identifying main flow pathways. *Computational Geosciences*, **23**(3), 617-629. https://doi.org/10.1007/s10596-019-9811-7

Strack, O. D. (2019). *Groundwater flow* (pp. 1-12). Encyclopedia of Water: Science, Technology, and Society.

Suk, H. (2017). Semi-analytical solution of land-derived solute transport under tidal fluctuation in a confined aquifer. *Journal of Hydrology*, **554**, 517-531. https://doi.org/10.1016/j.jhydrol.2017.09.033

Tang, D. H., Frind, E. O., & Sudicky, E. A. (1981). Contaminant transport in fractured porous media: Analytical solution for a single fracture. *Water Resources Research*, **17**(3), 555-564. https://doi.org/10.1029/wr017i003p00555

Toller, E. A. (2022). An analytic element model for intersecting and heterogeneous fractures in groundwater flow. *Water Resources Research*, **58**(5), e2021WR031520. https://doi.org/10.1029/2021wr031520

Tran, M., & Jha, B. (2020). Coupling between transport and geomechanics affects spreading and mixing during viscous fingering in deformable aquifers. *Advances in Water Resources*, **136**, 103485. https://doi.org/10.1016/j.advwatres.2019.103485

Turkman, M. A. A., Paulino, C. D., & Müller, P. (2019). *Computational bayesian statistics: An introduction* (Vol. 11). Cambridge University Press.

Vanderborght, J., Kasteel, R., Herbst, M., Javaux, M., Thiéry, D., Vanclillooster, M., et al. (2005). A set of analytical benchmarks to test numerical models of flow and transport in soils. *Vadose Zone Journal*, **4**(1), 206-221. https://doi.org/10.2136/vzj2005.0206

Vrugt, J. A., Gupta, H. V., Bastidas, L. A., Bouten, W., & Sorooshian, S. (2003). Effective and efficient algorithm for multiobjective optimization of hydrologic models. *Water Resources Research*, **39**(8). https://doi.org/10.1029/2002wr001746

Wang, W., Fan, D., Sheng, G., Chen, Z., & Su, Y. (2019). A review of analytical and semi-analytical fluid flow models for ultra-tight hydrocarbon reservoirs. *Fuel*, **256**, 115737. https://doi.org/10.1016/j.fuel.2019.115737

Wang, Y., Vuik, C., & Hajibeygi, H. (2022). CO2 storage in deep saline aquifers: Impacts of fractures on hydrodynamic trapping. *International Journal of Greenhouse Gas Control*, **113**, 103552. https://doi.org/10.1016/j.ijggc.2021.103552

Wei, X., Sun, J., & Sen, M. K. (2023). Quantifying uncertainty of salt body shapes recovered from gravity data using trans-dimensional Markov chain Monte Carlo sampling. *Geophysical Journal International*, **232**(3), 1957-1978. https://doi.org/10.1093/gji/ggac430
Werner, A. D., Bakker, M., Post, V. E., Vandenbohede, A., Lu, C., Ataie-Ashtiani, B., et al. (2013). Seawater intrusion processes, investigation and management: Recent advances and future challenges. *Advances in Water Resources*, **51**, 3–26. https://doi.org/10.1016/j.advwatres.2012.03.004

Witherspoon, P. A., Wang, J. S., Iwai, K., & Gale, J. E. (1980). Validity of cubic law for fluid flow in a deformable rock fracture. *Water Resources Research*, **16**(6), 1016–1024. https://doi.org/10.1029/wr016i006p01016

Wu, Y. S., Liu, H. H., & Bodvarsson, G. S. (2004). A triple-continuum approach for modeling flow and transport processes in fractured rock. *Journal of Contaminant Hydrology*, **73**(1–4), 145–179. https://doi.org/10.1016/j.jconhyd.2004.01.002

Xing, F., Masson, R., & Lopez, S. (2017). Parallel vertex approximate gradient discretization of hybrid dimensional Darcy flow and transport in discrete fracture networks. *Computational Geosciences*, **21**(4), 595–617. https://doi.org/10.1007/s10596-016-9606-z

Yih, C. S. (1961). Flow of a non-homogeneous fluid in a porous medium. *Journal of Fluid Mechanics*, **10**(1), 133–140. https://doi.org/10.1017/s0022112061000111

Yildirim, I. (2012). *Bayesian inference: Metropolis-hastings sampling*. Dept. of Brain and Cognitive Sciences, Univ. of Rochester.

Younes, A., Delay, F., Fajraoui, N., Fahs, M., & Mara, T. A. (2016). Global sensitivity analysis and Bayesian parameter inference for solute transport in porous media colonized by biofilms. *Journal of Contaminant Hydrology*, **191**, 1–18. https://doi.org/10.1016/j.jconhyd.2016.04.007

Younes, A., Hoteit, H., Helmig, R., & Fahs, M. (2022). A robust fully mixed finite element model for flow and transport in unsaturated fractured porous media. *Advances in Water Resources*, **166**, 104259. https://doi.org/10.1016/j.advwatres.2022.104259

Zaytoon, M., & Hamdan, M. H. (2021). Fluid mechanics at the interface between a variable viscosity fluid layer and a variable permeability porous medium. *WSEAS Transactions on Heat and Mass Transfer*, **16**, 159–169. https://doi.org/10.37394/232012.2021.16.19

Zhang, M., & Ayala, L. F. (2019). Variable rate and pressure integral solutions to the nonlinear gas diffusivity equation in unconventional systems. *Fuel*, **235**, 1100–1113. https://doi.org/10.1016/j.fuel.2018.08.065

Zhou, R., & Wang, Q. (2023). Analytical model for heat transfer in a discrete parallel fracture-rock matrix system. *Ground Water*, **61**(2), 183–192. https://doi.org/10.1111/gwat.13174

Zidane, A. (2023). CO2 viscosification for mobility alteration in improved oil recovery and CO2 sequestration. *Water*, **15**(9), 1730. https://doi.org/10.3390/w15091730