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

HAL
open science

Statistical physics and applied geosciences: some results
and perspectives

Benoît Noetinger

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

To cite this version:

Benoît Noetinger. Statistical physics and applied geosciences: some results and perspectives. Comptes Rendus. Physique, 2020, 21 (6), pp.539-560. 10.5802/crphys.40. hal-03135128v2

HAL Id: hal-03135128

https://ifp.hal.science/hal-03135128v2

Submitted on 25 Jul 2023

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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Comptes Rendus
Physique

Benoît Noetinger

**Statistical physics and applied geosciences: some results and perspectives**

Volume 21, issue 6 (2020), p. 539-560

Published online: 14 January 2021

Issue date: 14 January 2021

https://doi.org/10.5802/crphys.40

**Part of Special Issue:** Prizes of the French Academy of Sciences 2019 (continued)

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

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/

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

Les Comptes Rendus. Physique sont membres du
Centre Mersenne pour l'édition scientifique ouverte
www.centre-mersenne.org

e-ISSN : 1878-1535
<figure><img src="image_6.png" /></figure>

Prizes of the French Academy of Sciences 2019 (continued) / *Prix 2019 de l'Académie des sciences (suite)*

# Statistical physics and applied geosciences: some results and perspectives

## *Physique statistique et géosciences appliquées : quelques résultats et perspectives*

Benoît Noetinger<sup>a</sup>

**Abstract.** In this paper, some applications of statistical physics (SP) concepts and techniques in applied geosciences are reviewed. The domain includes hydrology, oil and gas industry, nuclear or CO2 waste geological disposal, massive energy storage, heat recovery from geothermal formations and many other applications. Several scales are considered: we start from applications of SP at the molecular scale to understand the effect of extreme confinements concerning the fluid transport in nanopores in clay rocks. The paper ends with coarse graining techniques that are employed to build operational models relevant at the practical scale of several kilometers, including strongly fractured geological environments. Perspectives are proposed regarding some issues about the practical use of these often over-parameterized models in connection to random matrix and graph theories and the associated quenched disorder problems and “big data” issues.

**Résumé.** Dans cette note, nous examinons quelques applications des concepts et outils de la physique statistique (PS) aux géosciences appliquées. Les applications vont de l'hydrologie au stockage de déchets nucléaires ou de CO2, au stockage massif d'énergie, en passant par l'industrie pétrolière ou gazière, et enfin les applications géothermiques. Vu la complexité intrinsèque des applications de terrain, les ingénieurs s'attachent en général à optimiser un critère économique, tout en veillant à maintenir la meilleure sécurité et en minimisant l'empreinte environnementale des projets. Il s'agit donc d'employer les connaissances les plus actuelles sur les transferts en milieu poreux. On s'intéresse à différentes échelles, des applications de la PS pour formuler les lois de transport dans des milieux extrêmement confinés tels des nanopores constituant la porosité des argiles. Ensuite, on présente les applications de la PS pour modéliser les écoulements à des échelles kilométriques intéressant les ingénieurs en charge d'application. On est dans une situation typique de désordre gelé, hors d'équilibre où les temps de relaxation peuvent être très longs, de l'ordre de plusieurs siècles. Des perspectives sont proposées pour utiliser des outils issus de la théorie des matrices aléatoires afin de faciliter l'utilisation pratique en "aide à la décision" de ces modèles bien souvent sur-paramétrés par des données elles-mêmes incertaines.

**Keywords.** Statistical physics, Applied geosciences, Porous media, Disorder, Nanopores, Upscaling, Quenched disorder.

**Mots-clés.** Physique statistique, Géosciences appliquées, Milieux poreux, Désordre, Nanopores, Changement d'échelle, Désordre gelé.

Available online 18th December 2020

<footnote><sup>a</sup> IFP Energies Nouvelles, France</footnote>

<footnote>E-mail: benoit.noetinger@ifpen.fr</footnote>
# 1. Introduction

In many applications of geosciences, the basic concern is to inject, to recover or to store some fluid, waste or even thermal energy in subsurface formations at good economical conditions, while minimizing overall environmental risks. As all these transport processes occur in the porosity of natural rocks, studying flow in porous media at several scales is a major fundamental issue. In order to do so, a concise macroscopic description of flow in porous media encompassing several scales without knowing the exact microstructure of the pore space must be obtained. That program is close in spirit to the program of SP, as it was acknowledged long ago by [1-4].

The paper will be focused on several applications of SP methods capable to treat subtle interface effects at the nanometer-scale, as well as understanding the effect of the quenched disorder of these porous rocks at large scales. It is well known that subsurface appears to be heterogeneous from nm to km scales [5]. That feature can lead engineers to build over parameterized numerical models that are used to help decision making on the field [6, 7]: well implantation, injection or recovery strategy, uncertainty management etc... Practical questions of major societal importance may be asked to anticipate the dispersion of some pollutant in an aquifer. How organizing a geothermal plant avoiding as far as possible induced seismicity having low social acceptance [8]? More long-term issues can be addressed: Global climate change may imply strong dryness and/or floods: how karstic formations will react to that change of solicitations is a question of major societal interest [9]. That can have dramatic consequences over water supply, and could imply catastrophic hazards. In that paper, we will illustrate some SP concepts that help to provide some methodology to answer such questions. These concepts arise essentially from the area of phase transitions and disordered systems, e.g. percolation theory. As the disorder is quenched, concepts and methods from spin glass theory may be useful. Self organization theories may also provide useful tools, especially for describing strongly non linear flows with retroaction between the flow and local rock transport properties. Finally, random matrix methods and random graphs are likely to provide useful concepts and results.

The present goal is to provide some answers to the following questions:

(i) is there some change of the analytical form of Darcy equation and of other transport equations if confinement effects due to small pore size arise? Are molecular simulations able to provide quantitative results, by correcting systematic bias due to the finite size of the simulation domain?

(ii) is there some change of the analytical form of the coarse grained Darcy equation and other transport equations as the support scale is changing?

(iii) Is there some self-averaging property and convergence to an effective homogeneous behavior of the system, and do we control this asymptotic convergence?

(iv) How to describe the large scale tracer advection/dispersion/diffusion in the associated velocity field, and more generally the mixing processes occuring in these random velocity fields?

(v) In the case of multiphase transport equations leading to hydrodynamic instabilitites such as Saffman-Taylor [10] (viscous fingering) or gravity driven instabilities, how does the development of the instabilities interact with the local conductivity heterogeneities?

(vi) A related issue concerns the strongly non-linear flow (such as non newtonian fluids, acidification leading to local changes of the conductivity, rock fracturation processes), can we observe some emergence of large scales patterns (fingers, wormholes) corresponding to some self-organised processes?

We discard in that paper any discussion about the application of SP to the evaluation of thermodynamics of bulk fluids required by applied geosciences as in many other applied science area. In addition, we do not discuss the increasingly popular Lattice-Boltzmann method that
originates from lattice gas automata (see [11] for the SP background of these methods) that is
used to upscale µCT scan images of rocks from pore to Darcy scale. A recent review may be found
in [12] and references therein. We will not discuss in depth the last item (vi), rich of essential SP
issues.

The paper is organised as follows: in next Section 1.1, we present some applications of molecular dynamics Monte Carlo studies that permit to get quantitative descriptions of flow in nanopores, even if important finite size bias must be accounted for properly. In a second section, we address the application of SP concepts to understand the form of averaged Darcy and transport equations in porous media characterized by random conductivity distributions Section 1.2. Averaging flow in fractured rocks that is an extreme case is described in Section 1.3 introducing percolation theory, dual porosity models that allow to describe these systems having highly contrasted relaxation times. Continuous time random walk techniques are presented, that yield very efficient computational techniques. In Section 1.4, we show how these considerations may be embedded using general random graphs and/or random matrix framework. Then some applications of random walks to tracer transport are given in Section 1.5. Then we give some elements regarding the influence of the disorder on two phase flows (Section 1.6) in which the disorder couples with the viscous fingering, yielding quite interesting ideas about the up-scaling of non-linear transport equations.

1.1. *Molecular dynamics at the pore scale*

Since the works of Darcy [13], developing a theory of flow through porous media from first principles was first viewed as a mathematical issue that can be treated by means of homogenization or volume averaging theories [14–17]. Such approaches are correct if the pore-scale Navier–Stokes description may be assumed on the bulk fluid. This description must be completed by the usual no slip boundary conditions at the surface of the pores, yielding a well posed problem. This description relies on the assumption that pore sizes are larger than a typical mean free path if gas transport is considered, or greater than several molecular diameters in the case of liquids. In that case, surface effects (with the notable exception of capillarity effects) may be neglected. In order to test the validity of this set of assumptions, a pioneering study was that of Koplik [18] in which Molecular dynamics (MD) techniques were reported to verify the validity of Navier Stokes description in a tiny pore. It was shown that the parabolic velocity profile as well as the Poiseuille relation relating the mean flow-rate to pressure drop remains quit robust, even for pore sizes of few molecular diameters. The longly debated question of the contact line motion between the rock, oil and water can be elucidated using a combination of MD and continuum methods [19]. The robustness of the usual non-slip boundary condition which combined with Stokes equation gives rise to Darcy’s law [20, 21] is more questionable. Experimental capacities, as well as improvements in SP description of non-equilibirum phenomena and continuous increasing power of computers, led to the so-called nanofluidics [22]. In present times, MD may be used to get a better understanding of moving contact lines too in realistic systems. The associated theoretical tools provided by SP allow then to propose rigorous coarse-graining procedures providing macroscopic law of interest for geoscientists [20, 21]. Note that a simplified version of MD gives rise to Lattice Boltzmann simulation algorithm in which molecular motions are over simplified, to give a fast and faithful description of the large scale flow. That allows to estimate directly quantitatively the permeability of rock samples using micro-scanner µCT images of the pore space [12].

Thanks to continuous increasing computing power, and to improvements in the characteriza-
tion of fluid and rocks, these MD techniques are presently adapted to model dynamics of com-
plex fluids confined between realistic mineral surfaces such as clays, the structure of which is
depicted in Figure 1.
<figure><img src="image_7.png" /><figcaption>**Figure 1.** Sketch of molecular dynamics simulations of transport between clay layers. The molecular structure of various clays is accounted for, as well as molecular motions under an imposed pressure gradient. Reprinted (adapted) with permission from [23]. Copyright 2020 American Chemical Society.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 2.** Diffusion coefficient of molecules in a confined pore of thickness *H*. On the right, the plain curve represents the result of the hydrodynamic calculation accounting for the systematic bias due to the infinite set of images of the unit cell. The dots are the set of values obtained by MD simulations for various *H*, showing the excellent agreement between simulation and theory. Reprinted (adapted) with permission from [23]. Copyright 2020 American Chemical Society.</figcaption></figure>

Such studies were at the origin of a better understanding of clay swelling phenomena [24, 25]. Density functional theories coupled with homogenization methods [26, 27] give quantitative descriptions of these complex natural media. Tiny electrostatic effects are amplified as soon as pore size may be compared with molecular dimensions, amplifying surface effects. These swelling effects are of primary importance to understand whether cap rock integrity is preserved in waste disposals or in CO₂ geological storage. If mechanical stresses due to swelling induce fracturation, the overall societal interest of the storage may be questionable. These successes motivated many SP studies to describe fluid motion at the nanoscale [20, 22, 25, 28–30]. The main goal is to be able to quantify the net effect of slippage and of electrostatic interactions at the pore boundaries, that may be neglected for usual pore sizes. But due to computing limitations,
these MD simulations are carried out on rather small unit cells usually supplemented by periodic boundary conditions. The convergence to large scale properties may be quite slow, inducing systematic errors. Analytical expressions are proposed to correct these systematic errors, both in the bulk or in nanopores, that is illustrated in Figure 2 [23, 31].

To conclude that section, SP tools such as MD or concepts (fluctuation theory, Kubo relations, density functional), coupled with continuous improvements of the knowledge of molecular fluid and rock properties are leading to rich applications that permit quantitative and more and more predictive descriptions [32]. This fine understanding of flow in nanopores and the associated changes in thermodynamical properties due to confinement have major consequences in applied geosciences, from nuclear waste disposal, CO₂ [33–35], to osmotic energy conversion using fresh and sea-water [36].

## 1.2. *From Darcy to large scale, up-scaling permeability and transport, anomalous dispersion*

We jump from nanopore scale to more macroscopic aquifer or reservoir kilometric scales to highlight applications of SP of disordered systems to geosciences, even if SP concepts and methods such as percolation theory were usefully employed at the scale of the cores (some cm). In Figure 3, we sketch the overall scales of reservoir (or aquifer) numerical simulations by depicting the characteristic sizes that enter in any aquifer description in general followed by a numerical simulation because very few analytical solutions are available.

Most applications share the following issue: solving a diffusion-like equation that reads:

$$ \phi c_t \frac{\partial p(\mathbf{r}, t)}{\partial t} = \nabla \cdot (k(\mathbf{r}) \nabla p(\mathbf{r}, t)) + f(\mathbf{r}). \quad (1) $$

Here, the parameters $k(\mathbf{r})$, $p(\mathbf{r}, t)$ and $f(\mathbf{r})$ denote respectively the local conductivity, time-dependent potential and a source term. Dirichlet or Neumann Boundary conditions are known at the boundary of the domain. All these quantities depend on position vector $\mathbf{r}$ and time $t$. In most cases, the quenched positive $k(\mathbf{r})$, is represented as being a random function of position. It is characterized by some mean-value and fluctuations that are measured by the geologist, by means of the so-called geostatistical approach, originally founded by Matheron using statistical concepts [37–39]. Log-normal distributions (the logarithm of the conductivity is a Gaussian distributed variable) were observed at the core scale in well defined geological environments [40]. Note that in most cases, the amount of data is not sufficient to provide the form of the probability distribution function (pdf) of $k$, and high order correlation functions are impossible to determine. As a consequence, even the input stochastic model is questionable. This explains the quasi infinite set of approaches that exist to generate random fields compatible with geological observations, the textbook [41] present the most popular approaches and methods.

The basic issue is to be able to quantify the net effect of heterogeneities about the behaviour of the aquifer, as well as the related uncertainties. In that context, such issue were investigated by hydrogeologists [40, 42–44], mathematicians [45, 46] and physicists, among other [47–54]. So the basic issue is to transfer the small scale spatial fluctuations to a large scale support that encompasses the low spatial frequency components of the fields of interest. The calculations are carried out in practice using some numerical model in which the Laplace equation is solved using a grid of resolution $L$ generally considerably much coarser than the input fine geological grid $\Delta$ (notations in Figure 3), because the available computing power leads also to continuous improvement of local geological 3D representations. This implies obtaining a coarse grained Laplace equation with a renormalized conductivity map accounting as best as possible to the local subgrid variations. This is a classical issue addressed long ago by Maxwell, Landau and Lifshitz among others, [40, 55–57]. There is a great amount of literature using SP concepts
such as percolation theory, real-space renormalization techniques, field theoretical methods including diagrams summation techniques [41, 47–49, 53, 54, 58–61]. A difficulty is to merge these sophisticated techniques with the pragmatical needs of field engineers solving real time issues. One can average the solution of (1) to get the average head (or pressure) $\langle p(\mathbf{r}, t) \rangle$ in which the ensemble-average $\langle \cdots \rangle$ is to be taken on the quenched disorder of the conductivity field $k(\mathbf{r})$.

Once the velocity field is given, advection dispersion diffusion equation may be solved, e.g. to anticipate the motion of some pollutant, or heat recovery. Modelling multiphase flow leads to solve non linear equations involving the repeated solving of Laplace equation (1), as described in Section 1.6.

### 1.2.1. *Averaging Darcy's law*

It can be shown that under quite general hypothesis (statistical stationarity and convergence conditions) that the average potential $\langle p(\mathbf{r}, t) \rangle$ is driven by an effective equation that reads [61]:

$$ \phi c_t \frac{\partial \langle p(\mathbf{r}, t) \rangle}{\partial t} = \int_{-\infty}^{t} dt' \int d^D \mathbf{r}' \nabla \cdot (\Sigma(\mathbf{r} - \mathbf{r}', t - t') \nabla p(\mathbf{r}', t')) + f(\mathbf{r}). \quad (2) $$

The average local flux is a spatially weighted time convolution of the average potential gradient around the considered time and location. It means that the underlying disorder couples different points. That is reminiscent of the “overlap” that arises in quenched average that are provided by replica methods in spin glasses theories, that are the archetypel of quenched disorder approaches in SP [62]. The kernel $\Sigma(\mathbf{r}, t)$ may be obtained as a result of a summation of 1P irreducible irreducible diagrams of a perturbation expansion of the solution of (1), in a power series of the conductivity fluctuations, a so-called self-energy [61]. The diagram resummation techniques familiar to SP allows to build in a systematic manner from the perturbation expansion of the average solution $\langle p(\mathbf{r}, t) \rangle$ the equation driving $\langle p(\mathbf{r}, t) \rangle$. As the Green's function of Laplace operator is long-ranged, it is far more practical to manipulate an effective equation. The kernel $\Sigma(\mathbf{r}, t)$ involves a rather complex series of integrals involving the heat kernel and correlation functions of the conductivity of higher and higher order. In the generic case, the spatial range of $\Sigma(\mathbf{r}, t)$ is controlled by the correlation length of the conductivity fluctuations, and it time range is a typical diffusion time over this correlation scale. This structure explains the self-averaging character of the Laplace equation: low frequency components of the potential on one single large realization of the domain behaves like a Monte Carlo average of the potential over many independent realizations of the disorder [46, 61]. This explains why pumping tests (corresponding to point-wise solutions) “homogenize” by themselves It may be shown that for an infinite domain, at long times, the average potential is driven by the following equation:

$$ \phi c_t \frac{\partial p(\mathbf{r}, t)}{\partial t} = \nabla \cdot (K_{\text{eff}} \nabla p(\mathbf{r}, t)) + f(\mathbf{r}). \quad (3) $$

The parameter $K_{\text{eff}}$ given formally by $K_{\text{eff}} = \int_{0}^{+\infty} dt \int d^D \mathbf{r} \Sigma(\mathbf{r}, t)$ is called the effective conductivity. It corresponds to the “natural” large scale relation between the mean flux and the large scale pressure gradient that can be provided by homogenization theories by means of the so-called “auxiliary problem” to be solved numerically in x, y and z directions [46] and references therein, as it is illustrated Figure 4. It can be shown that at long times, for any realisation of the disorder, the behavior of the system will converge to this equation: this is another manifestation of the self averaging property [46].

Many investigators proposed expressions relating $K_{\text{eff}}$ to the underlying disorder [64]. For 1 dimension, an elementary analytical calculation provides the harmonic average $K_{\text{eff}} = \langle k^{-1} \rangle^{-1}$. For $D=2$, a nice duality argument [14] shows that $K_{\text{eff}} = \exp\langle \text{Log}(k) \rangle$ if the conductivity distribution is log-normal. No simple general analytical expression exists in the general case for $D \ge 3$, even if
<figure><img src="image_9.png" /><figcaption>**Figure 3.** Geometry of the problem on a simplified 2D section of an aquifer or reservoir model domain $\omega$. A working coarse grid made of blocks $\Omega$ of typical size $L$ is superimposed to a geological fine grid of typical size $\Delta$ in order to solve the mass conservation equations of interest. The intermediate scale $\lambda$ serves in posterior treatments for checking the overall consistency of the model. Reprinted (adapted) from Ref. [63]. Copyright 2020 with permission of Elsevier.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 4.** Up-scaling geometry. The coarse block of size $L$ have a detailed conductivity map given by the geologist. It is up-scaled by solving a steady-state quasi Laplace equation to determine an effective conductivity in the mean flow direction that will serve as input of the simulator at coarse scales. Changing the direction of the mean driving pressure gradient allows to determine a conductivity tensor. The boundary conditions are usually no-flux parallel to the imposed mean flow, or periodic.</figcaption></figure>

a great deal of research was devoted to develop such analytical expressions using additional hypothesis. SP techniques such as field theoretical methods [47, 48, 53, 60, 63] were employed to find some robust approximations. Many authors attempted to justify the so called Landau–Lifschitz–Matheron (LLM) [14, 56] formula that reads:

$$K_{\text{eff}} \approx \langle k^{(1-\frac{2}{D})} \rangle^{\frac{1}{(1-\frac{2}{D})}}. \quad (4)$$
<figure><img src="image_11.png" /></figure>

**Figure 5.** Monte Carlo study of the evolution of the effective conductivity pdf with coarsening scale $\lambda$. The overall flow is solved using several realizations of the input log conductivity map (left). The associated local dissipation map (center) allows to evaluate a distribution of coarsened effective conductivities averaged at scale $\lambda$ [63]. The resulting pdf's are plotted (right). The self averaging (homogenization) is highlighted by the sharply peaked distribution around the geometric average for $\lambda = 128$ units, the stability across scales of the Gaussian distribution may be observed too. Reprinted (adapted) from Ref. [63]. Copyright 2020 with permission of Elsevier.

This formula is found to be exact in 1 and 2D and up to fourth-order in the log-conductivity variance [65]. It remains an approximation in 3D or more [43, 54, 66], moreover, it is quite robust when compared with numerical tests in the case of log normally distributed conductivities [49, 63, 67]. Series resumption techniques [48], renormalization group (RG) arguments [48–50] give some clues in favor of this robustness. The author still thinks that there is some hidden powerful theoretical framework to be developed justifying its robustness and giving some sense to this formula. In 2 dimensions, numerically the log normal distribution appears to be stable under the up-scaling transformation, Figure 5 from [63], analogous to the central limit theorem. That may be related to the duality argument of Matheron [14] that justifies the geometric mean in 2D. In the 3D case, studying the emergence of a “stable” conductivity distribution invariant on the up scaling (RG) transformation would also be useful for studying strongly correlated systems having conductivity correlations decaying as a power law with the lag distance. A related question is to give some sense to the so-called uncorrelated case which may only be valid at a given observation scale that plays the role of the fixed scale while letting the ultraviolet cut off going at infinity, but keeping observed quantities fixed to their nominal values [68]. This may be illustrated on Figure 5 below: in the practical side, in most cases, at present times, using a numerical technique is sufficient. In case of extremely heterogeneous media (that can correspond to bimodal media) at percolation threshold, the well-known percolation transition may occur [58, 59, 63, 69, 70]. The percolation second order transition may be observed on Figure 6.

## 1.3. Fractured rocks

### 1.3.1. Upscaling the fracture network

Rocks are almost always fractured, and it may be expected that these fractures control the flow. In this paper, self-organization effects arising from geomechanics (a SP issue in itself!) controlling the overall organization of fracture networks will be ignored [71]. Understanding such flows is essential in many applications, such as geothermal applications, water resources management and oil and gas recovery. That explains the many *ad-hoc* rather empirical models that were developed, in particular the popular double porosity model that represent the rock and the fractures as two superimposed continua [72, 73].
<figure><img src="image_12.png" /><figcaption>**Figure 6.** Monte Carlo study of the evolution of the effective conductivity pdf with coarsening scale $\lambda$. The overall flow is solved using several realizations of the input bimodal map (left). The associated local dissipation map (center) allows to evaluate a distribution of coarsened effective conductivities averaged at scale $\lambda$ [63]. The resulting pdf's are plotted (right). As the scale is increasing, the two peaks merge, the bimodal distribution disappears and becomes a (log-normal like?) distribution. The convergence to this asymptotic distribution shows critical slowing-down when the facies proportions are close to percolation threshold. In the infinite contrast case, scaling-laws are recovered [70]. Reprinted (adapted) from Ref. [63]. Copyright 2020 with permission of Elsevier.</figcaption></figure>

In the case of fractured rocks, the heterogeneity is extreme: conductivity may vary by several orders of magnitude between the matrix (that stores the quantity of interest, fluids, energy etc...) and the fractures of almost vanishing measure which carry most of the flow to the outlet. These fractures may have random orientations, power-law distributions of lengths an apertures [74–76] leading to quite complex parameter space and phase diagrams. In that situation, continuous perturbation theories break-down, and other methods must be employed. Percolation theory approaches focused on the role of the fracture network connectivity [58, 59, 70, 74, 76–79] provide an excellent framework for describing these systems controlled by connectivity effects, at least if the system is close to the percolation threshold.

Fracture networks can be quite naturally represented as a random resistor network which corresponding graph shares vertices that represent the fracture intersections and edges representing the connection (fractures) between these intersections. Although the mapping is straightforward in 2D [66, 80], it is far more complex to derive it rigorously in the 3D case considering fractures that are 2 dimensional objects [81, 82] embedded in the usual 3D space. The resulting random resistor network corresponds to a low order approximation that can be improved systematically. This random resistor network may be associated to a weighted Random graph. The associated Laplacian matrix summarizes the hydraulic connections between fractures. Random graphs and random matrices theoretical techniques could provide a useful framework [83–93].

### 1.3.2. Coupling the fracture network with the matrix, beyond the dual porosity model

Coupling the fracture network with the matrix characterized by larger relaxation times is generally done using the classical dual porosity model of Barenblatt et al. [72]. This model allows to account for the smallest relaxation time of the matrix that is related to the eigenvalues of the Laplace operator acting on the matrix domain with Dirichlet boundary conditions [94]. More refined models attempted to improve that description by adding several relaxation times [94–96]. Continuous time random walk (CTRW) techniques were proposed for fast computing of the relaxation of the matrix, thanks to a direct relation between the first exit time distribution of a particle undergoing Brownian motion in the matrix, and the matrix relaxation function directly related to the residence time distribution in the matrix [82, 97, 98]. The same method allows to
<figure><img src="image_13.png" /><figcaption>**Figure 7.** A percolation network $p = 0.57$, red = isolated clusters, black = backbone, green = dead-ends. The zoom on the right shows the underlying matrix in which CTRW are performed, in white. Reprinted (adapted) from Ref. [99]. Copyright 2020 with permission of Springer-Nature.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 8.** Exchange coefficient computed for the percolation network $p = 0.57$. Reprinted (adapted) from Ref. [99]. Copyright 2020 with permission of Springer-Nature.</figcaption></figure>

compute the effective conductivity using Einstein relation for the mean square displacement of the diffusing particle inside the fracture domain. Such a method was implemented on fracture networks generated as 2D bond percolation Figure 7 [99].

In Figure 8, we plotted the dependence $p - p_c$ of the mean residence time $\langle t \rangle$ in the matrix (of diffusivity one unit) with respect to the proportion of active fractures. The different set of points correspond to different mean residence times $\langle t_{exit} \rangle$ inside the matrix that depend on the fracture subnetwork that was kept, i.e., all the fractures including non-relevant isolated clusters, or only the percolation backbone without dead ends. Intermediate curves correspond to different
treatments of the remaining clusters. One can note that keeping the whole set of fractures (red dots) (Figure 8) does not lead to any critical divergence of the mean residence times close to $p_c$. Further studies must be carried out to get a better characterization of the associated critical exponents. It can be remarked that these mean residence time that can have different values may represent relaxation associated to different transport models (pure diffusion in the matrix, diffusion in the fractures, or advection diffusion in the fractures). In practice, it means that the coupling coefficient determined by pressure tests (pure diffusion) on a given fractured reservoir may not be directly applicable to geothermal applications (advection dispersion in the fractures, heat diffusion in the matrix), that concerns the fracture network backbone.

## 1.4. Up-scaling, graph and random matrix theory

Fracture networks appear naturally as random graphs that provide Laplacian matrices In 2D cases, in which fractures may be viewed as 1D tubes, solving the dominant flow in the fracture network is strictly analogous to determine currents in a random resistor network, the nodes of which being the intersections between fractures [80]. It can be shown that such a construction can be carried out in 3D, even if the intersections between fractures are segments [81, 82]. So, one is led to solve large linear systems of equation that reads

$$ \forall i \sum_{j \in \langle i \rangle} T_{ij} (P_j - P_i) = Q_i \qquad (5) $$

$$ \sum_{i} P_{i} = 0. \qquad (6) $$

The source terms $Q_i$ are such that $\sum_i Q_j = 0$. The set of labels $\langle i \rangle$ denotes the set of vertices connected to vertex $i$ by one edge. Note that it is also the case for discretized equations corresponding to Darcy flow in heterogeneous random systems discussed in Section 1.2. Considering the discrete equation corresponding to a Darcy problem, the set of $\langle i \rangle$ is essentially the 2D neighbours of a given cell (using other numerical techniques will essentially change this set of neighbours and the values of the $T_{ij}$). The coupling coefficients $T_{ij}$ are related to the underlying conductivity maps [63]. In the fractured case, the set of neighbours $\langle i \rangle$ may be arbitrary large, so in that situation, two superimposed disorders are superimposed, one from the structure of the graph, the other from the $T_{ij}$'s. Gathering all the unknowns in N dimension vectors, the equations may be written under a more compact form $\mathbf{A} \cdot \mathbf{P} = \mathbf{Q}$. The operator $\mathbf{A}$ appears as a random matrix corresponding to a weighted graph of an associated random conductivity network [85, 86, 100–103]. In the fractured case, these matrices may be treated using methods of graph theory [92, 93, 104, 105].

It appears interesting to study the distribution of the “small” eigenvalues of $\mathbf{A}$, (as well as the corresponding eigenvectors), the null value having a multiplicity equal to the number of connected components of the associated graph [83]. Retaining one component, these eigenvalues may be denoted by $\lambda_1 = 0 \le \lambda_2 \le \lambda_q \le \lambda_{N_c-1}$ n which $N_c$ denotes the number of vertices of the retained connected component. In order to illustrate the idea, consider the case of a simple path graph (node $i$ connected to nodes $i-1$ and $i+1$, at the exception of nodes 1 and $N_c$ having only one connection). The corresponding linear system (6) can be solved easily by recursion, its solution may thus be averaged over the disorder of $T_{ij}$. This solution is itself given by the solution of an effective linear system sharing the same structure than the original one equation (6) by setting $T_{eff} = \langle T_{ij}^{-1} \rangle^{-1} = \langle T^{-1} \rangle^{-1}$. So in that simple case, the effective set of equations $\mathbf{A}_{eff} = \langle \mathbf{A}^{-1} \rangle^{-1}$ has a spectrum given by $\langle T^{-1} \rangle^{-1} 4 \sin^2(\pi q / 2N_c)$. That spectrum is itself equivalent to $[\pi^2 / N_c^2] q^2$ (small $q$) [83]. The reader should note that $\mathbf{A}^{-1}$ has a well defined sense working on properly defined subspaces. In real space domain problems, small $q$ corresponds to small frequencies (large lengthscales). But studying small eigenvalues keeps a well-defined sense without being embedded in an Euclidean framework. On that aspect, homogenization theory states that for small $q$, a Laplacian
operator with oscillating coefficients behaves may be replaced by an effective laplacian with an effective conductivity $T_{\text{eff}}$ [46], with eigenvalues scaling as $T_{\text{eff}}[4\pi^2/N^2]q^2$ (small $q$). Studying the small eigenvalue spectrum combined with self-averaging properties may lead to new results concerning up scaling: low frequency eigenvalues may behave as $T_{\text{eff}}q^2$. This implies that in the case of a random graph with random weights $T_{ij}$, by identification of the distribution of the smallest eigenvalues, it could be possible to get an effective conductivity $T_{\text{eff}}$, and an effective dimension $D$ (possibly > 3 in case of highly connected media) that will characterize the average connectivity structure of the problem at hand. In particular, is it possible to propose a simple averaging formula such as $T_{\text{eff}} \approx \langle T^{(1-2/D)} \rangle 1/(1-2/D)$? In the considered case, the so-called uncorrelated conductivity distribution has a well-defined sense, in opposition to the corresponding continuous limit in which this concept is basically meaningless [63, 106]. Assuming uncorrelated $T_{ij}$ sharing the same probability distribution, $A_{\text{eff}}$ must be related to the adjacency matrix of the graph of the average Laplacian matrix $\langle A \rangle$. Considering the associated resolvent and Stieltjes transform will provide information about all the moments of the random matrix and also about its density of eigenvalues in the large $N$ limit. In particular, the spectrum of such large matrices that depends randomly on the geological input parameters may provide information about the effective number of relevant degrees of freedom that must be retained to describe the system [107]. Depending on the degree of disorder, some eigenvalues of the associated Laplacian matrix can be expected to provide information while other eigenvalues may follow some universal distribution such as a Marchenko-Pastur distribution [108]. Such topics are deeply connected to classification methods by neural nets [109, 110].

## 1.5. Transport and mixing issues for passive and reactive flows

A classical issue is to be able to model the spreading of a passive tracer undergoing advection and diffusion/dispersion in a given steady state imposed flow in an heterogeneous/fractured medium. That issue motivated many approaches since the early approach of Saffman [3]. The porous medium was represented as a network of tubes connected at several pore intersections. The tracer spreading was thus represented as deterministic motion of the tracer particles inside the pores, being randomized by random choice of the pore at the intersection. This gives rise to a macroscopic dispersion, the overall motion of a cloud of tracer being described by the following equation:

$$ \phi \frac{\partial C(\mathbf{r}, t)}{\partial t} + \nabla \cdot (\mathbf{U}(\mathbf{r}, t) C(\mathbf{r}, t)) = \nabla \cdot (\mathbf{D} \nabla C(\mathbf{r}, t)). \quad (7) $$

In that equation, $\mathbf{U}(\mathbf{r}, t)$ denotes the Darcy velocity, which is divergence free, implying strong range spatial correlations. A great deal of efforts were carried-out at pore scale, including SP techniques, in order to relate the value of tensor $\mathbf{D}$ to some descriptors of the velocity field and to the molecular diffusion coefficient [111, 112]. The basic mechanisms include Taylor dispersion [113, 114], and the amplification of spreading due to the presence of stagnation point of the flow. These points play the role of bifurcation like zones that over amplify local processes such as smaller scale diffusion/dispersion can be pointed out [115–117]. Specific random walk simulation techniques allow to estimate breakthrough times (BTC) for various Péclet numbers [99, 117], as illustrated in Figure 9 below: at a larger scale, stochastic hydrology [40, 42] techniques were developed to obtain an up-scaled description of the tracer motion. The approach starts from a Darcy scale description that includes an input dispersion tensor $\mathbf{D}_{\text{eff}} \approx \langle U \rangle \sigma_{\text{log},k}^2 l_c$ that accounts for subscale effects, that is the so called macrodispersion phenomenon that yields a large scale Fickian like description, as soon as the spatial correlations of the conductivity are short-ranged. In that context, the second order moment of the tracer spreading grows linearly with time. In another study, Matheron and de Marsily [118] shown that transverse diffusion effects on a stratified
<figure><img src="image_15.png" /><figcaption>**Figure 9.** Left: three-dimensional visualization of the pore network (blue) and the connected microporosity (yellow) for a sub-volume of 1203 voxels of the (0.8 mm³) Berea sandstone sample studied by [117]. Right: BTCs for pore-scale mobile-immobile transport through the Berea sample for different values of the Peclet number computed by TDRW, Reprinted (adapted) from Ref. [99]. Copyright 2020 with permission of Springer-Nature.</figcaption></figure>

flow may lead to non-diffusive over all behaviour for the longitudinal dispersion. This surprising phenomenon is mainly due to the strong probability of return to the origin regarding transverse diffusion that creates long time tail correlations leading to superdiffusion in the direction of stratification that has an infinite correlation scale. Many contributions, including the one of [119, 120] led to more general descriptions. In another study, renormalization-group ideas “à la Wilson” were set-up to compute the spreading of a tracer in a heterogeneous medium [121]. Field applications provide encouraging results. At the pore and core scale, anomalous dispersion effect leading to a description in terms of fractional derivatives were proposed [122–126]. NMR techniques allow to determine the “anomalous” parameters from laboratory tracer experiments. In the case of fractured media, approaches combining anomalous transport concepts and multiple media ideas may be combined to advanced continuous time random walk techniques [92, 93, 95, 99, 127]. This results in the possible emergence of robust stable laws that can be tested on real field experiments. A major related issue of interest that must be addressed in the case of reactive flows is to describe the intimate mixing occurring between fluids flowing in porous media, in which the useful concept of lamella diffusion appears to be very promising [27, 128–131].

## 1.6. *Non-linear issues, coupling between the quenched disorder and flow instabilities*

In this section, we are interested by strongly non-linear processes that arise in two phase flow displacements by fluids having different mobilities in the rock. This can be the case of water displacing oil, of CO₂ injection etc... This issue gives rise to the well-known Saffman–Taylor instability called viscous fingering if the displacing fluid is more mobile than the fluid initially in place [132]. Due to its importance for the secondary oil recovery applications (displacing oil with water), this seminal work gave rise to many subsequent works [10, 133–135]. Many SP concepts such as diffusion limited aggregation (DLA) corresponding to the extreme case of air displacing a liquid, fractals were illustrated by experiments involving fingering [10]. The question of the selection of the ultimate finger pattern received many attention [136] in the eighties.
Many of these studies were performed in Hele-Shaw cells on in well-controlled micromodels. Analogous phenomena occur in natural porous media, and the stability criteria can be adapted using the standard description of two phase flow using relative permeability concepts in 3D rocks, as it was explained in the relatively ignored paper of King and Dunayevsky [137] and references therein. It is well known that in that case, it appears a well defined sharp front separating the fluid originally in place and a rich phase that moves at the local velocity of the fluid, as shown in Figure 10, top. That front may become in turn unstable, by the same mechanism than the Saffman–Taylor instability, as shown in Figure 10, bottom. That phenomena arises from the coupling between the fluid motion and the up dating of the local mobility, the criterion being the so-called total mobility jump evaluated ahead and behind the front. The finger selection process is mainly controlled by the underlying heterogeneity of the rock. It is easy to understand that a highly mobile fluid will follow high velocity paths, its presence amplifying thus that advantage by a positive feed-back loop. This is the so called channeling issue. This problem was addressed by De Wit and Homsy [138, 139], and revisited in the stochastic context by [140–142]. The idea of the latter contributions was to adapt the theory developed by King and Dunayevsky [137] using single phase flow perturbation theory techniques. The underlying equations reads:

$$ \nabla \cdot [\lambda(S(\mathbf{r}, t))k(\mathbf{r})\nabla p(\mathbf{r}, t)] = 0 \qquad (8) $$

$$ \phi \frac{\partial S(\mathbf{r}, t)}{\partial t} + \nabla \cdot (f(S(\mathbf{r}, t))\mathbf{U}(\mathbf{r}, t)) = 0 \qquad (9) $$

$$ \mathbf{U}(\mathbf{r}, t) = -\lambda(S(\mathbf{r}, t))k(\mathbf{r})\nabla p(\mathbf{r}, t). \qquad (10) $$

Here, $(S(\mathbf{r}, t))$, $\lambda(S(\mathbf{r}, t))$ and $f(S(\mathbf{r}, t))S$ denote respectively the water saturation (local % of water, total mobility and the so-called fractional flow of water). This set of coupled equations may be solved numerically (it is at the heart of any multiphase flow in porous media simulator, to which additional complexities such as phase transitions and boundary condition management must be added). The saturation equation (9) is hyperbolic, leading to the formation of a shock front, whose stability is controlled by the jump of the total mobility $\lambda(S(\mathbf{r}, t))$ at the front.

The technical difficulty for setting-up a perturbation expansion comes from the presence of the front that implies a mobility jump that renders perturbation theory a bit tricky [137, 141]. The difficulty may be avoided by a suitable change of variable, using a working variable $x(S, y, t)$ rather than $S(x, y, t)$. It is thus possible to introduce the function $x(S_f, t)$, in which $S_f$ is the saturation of water just behind the front $M_f$ is the corresponding total mobility jump $M_f = \lambda(S_f)/\lambda(S = 0)$. The randomness of the underlying conductivity field propagates to the randomness of $x(S_f, t)$, of average value $\langle \mathbf{U} \rangle_t$. At long times, the associated two point correlation function can be shown to converge to a well defined function in the stable case, while in unstable case it diverges, a manifestation of the spreading of the front (even if some logarithmic singularities are remaining, due to the singular character of the instability at large wavelengths [143, 144]). A possible approach of practical interest close to the single phase flow approach would be to look at an effective equation driving the ensemble-averaged water saturation $\langle S(\mathbf{r}, t) \rangle$, or the Y-averaged saturation $S(x, t)$. A diffusive regime arises in the case $M_f = 1$, that corresponds to a marginal stability criterion (corresponding to an order parameter of a phase transition), leading to a macrodispersion equation similar to (7). In the general case, several proposals were reported long ago for characterizing the emerging large scale transport equation [145–149]. In the unstable case, one can consider that long fingers parallel to the imposed flow may be treated as a stratified medium. This leads to modify the fractional flow function with an ad-hoc change [145]. In the stable case, it can be shown that the competition between the disorder that distorts the front and the viscous forces that tends to sharpen the front [141] must lead to another form of the effective fractional flow, including some macrodispersion representing the net effect of the averaged disorder. In the infinite
<figure><img src="image_16.png" /><figcaption>**Figure 10.** Simulated dynamics of a two phase flow front in a heterogeneous rock, imposed mean flow from left to right, underlying random log normally distributed conductivity map top (a) stable case, (b) Y averaged saturation at different times, bottom (c) fingering in the unstable case, (d) associated Y-averaged saturation at different times.</figcaption></figure>

contrast case (e.g. immiscible gas injection) diffusion limited aggregation (DLA) models were pro-
posed, leading to a very rich literature involving percolation invasion models, fractals [150–153]
with many contributions of SP.

1.7. *Conclusions and some perspectives*

In that short review, some connections between statistical physics approaches and applied geosciences were presented. These connections are not new and are continuously enriching both communities. At the molecular scales, SP techniques including molecular dynamics tools help to find the form of the constitutive relation relating fluxes to gradient of potential, even for charged real fluids in extremely confined environments. At pore scale, SP may help to describe complex multiphase flow, in particular [19] and by improving lattice Boltzmann methods.

In practice, an essential feature of these natural systems is their overall insufficient charac-
terization, implying that the modeller must find the optimal balance between the details of the
model, and the lack of information. A very beautiful model, but completely over-parameterized
can be useless [73]. So, one of the first task of the modeller is to be able to identify flow regimes,
and the aggregates of most relevant input parameters by means of phase diagrams highlight-
ing the most relevant parameters that control the overall behaviour of the system. Geoscientists
must manage huge parameter space and the propagation of the uncertainties due to the lack of
an exhaustive description of the system. SP tools help to understand the role of the quenched
disorder of the medium present at all scales. Drawing the “phase diagram” of the problem at
hand, i.e. the set of dominant parameters controlling the overall behaviour of the system, and
providing descriptions of the critical behaviour of the system between the different regions of the phase diagram. Providing coarse renormalized equations describing averaged potential spatial variations, or tracer spreading moments with accuracy with few parameters is a promising approach [154, 155]. Studying the strong coupling (competition or amplification) between the disorder of the porous medium, and the development of viscous and gravitational instabilities such as Saffman–Taylor viscous fingering remains a deep issue. Focusing on the front dynamics moving in a random medium using some KPZ like approaches [156] could be a rich research avenue. A combination of the methods developed by King and Dunayevsky [137, 141] and of stochastic perturbation theory could be an interesting approach, although the existence of a local equation driving the front dynamics remains questionable, in view of the long range character of the Laplace equation Green's function. Such approaches could help to improve empirical descriptions [132, 141, 142, 149, 156, 157].

Studying self-organization phenomena that arise once the transport processes at hand modify the porous media structure with a strong feed-back, implying a strong non-linear coupling [71, 158, 159] using SP concepts is a major research avenue. For engineers, this may help to constrain the parameter space of the problem at hand, in particular the stochastic models of discrete fracture networks that may be over parameterized. These models are close to morphogenesis models of SP [160–162] that may help to describe overall fracture network organization. Such models may be relevant in the context of geomorphology, that could in turn provide some information about the statistical properties of the quenched disorder that was discussed throughout the paper. Those topics involving advanced geoscience and SP concepts, coupled to global climate evolution are well beyond the scope of present paper [153, 163–165].

Subsurface is intrinsically a quenched system that falls in the area of spin glasses issues [62]. The considerable number of degrees of freedom made it a natural candidate for using big data and/or artificial intelligence techniques in order to help the engineers to manage their intrinsic complexity, and to select the most relevant approaches and parameters. As it was suggested in Section 1.4, studying the spectrum of such large random matrices depending on the geological input parameters may provide information about the effective number of relevant degrees of freedom that must be retained to describe the system [107]. Depending on the degree of disorder, some eigenvalues of the associated Laplacian matrix can be expected to provide information while other eigenvalues may follow some universal distribution such as a Marchenko–Pastur distribution [108]. Such topics are deeply connected to classification methods by neural nets [109, 110]. This set of techniques may be relevant for solving the inverse problems (modifying the model parameters or the model itself to account for continuously arriving data). Such inverse problems may be solve by minimizing a suitable error functional accounting for both these data and prior information [166–169]. All these topics are connected to each other by means of recent developments of spin glass theories [62, 170] that provide useful approaches and algorithms to provide rigorous and operational foundations to uncertainty management allowing to make the best decision with large random parameter sets associated with complex physics. We expect that this short overview highlights numerous applications of statistical physics to applied geosciences and that it will stimulate discussions to build bridges between very active areas of research in statistical physics and geosciences.

## Acknowledgments

The author is indebted to Mrs Mickaële le Ravalec, Pauline Simonnin and MM Vincent Artus, Alejandro Boschan, Ivan Colecchio, Fred Delay, Marco Dentz, André Fourno, Philippe Gouze, Tri Dat Ngo, Alejandro Otero, Ludovic Ricard, Benjamin Rotenberg for their useful contributions.
References

[1] A. E. Scheidegger, "Statistical hydrodynamics in porous media", *J. Appl. Phys.* **25** (1954), no. 8, p. 994-1001.

[2] A. E. Scheidegger, "The physics of flow through porous media", *Soil Sci.* **86** (1958), no. 6, p. 355.

[3] P. Saffman, "A theory of dispersion in a porous medium", *J. Fluid Mech.* **6** (1959), no. 3, p. 321-349.

[4] A. E. Scheidegger, E. F. Johnson, "The statistical behavior of instabilities in displacement processes in porous media", *Canad. J. Phys.* **39** (1961), no. 2, p. 326-334.

[5] A. Guadagnini, M. Riva, S. P. Neuman, "Recent advances in scalable non-Gaussian geostatistics: the generalized sub-Gaussian model", *J. Hydrol.* **562** (2018), p. 685-691.

[6] S. Gorell, R. Bassett *et al.*, "Trends in reservoir simulation: big models, scalable models? Will you please make up your mind?", in *SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers*, 2001.

[7] F. J. Floris, M. Bush, M. Cuypers, F. Roggero, A. R. Syversveen, "Methods for quantifying the uncertainty of production forecasts: a comparative study", *Petrol. Geosci.* **7** (2001), no. S, p. S87-S96.

[8] N. Deichmann, D. Giardini, "Earthquakes induced by the stimulation of an enhanced geothermal system below Basel (Switzerland)", *Seismol. Res. Lett.* **80** (2009), no. 5, p. 784-798.

[9] H. Loáiciga, D. R. Maidment, J. B. Valdes, "Climate-change impacts in a regional karst aquifer, Texas, USA", *J. Hydrol.* **227** (2000), no. 1-4, p. 173-194.

[10] P. Tabeling, G. Zocchi, A. Libchaber, "An experimental study of the Saffman-Taylor instability", in *Dynamics of Curved Fronts*, Academic Press, 1988, p. 219-234.

[11] D. H. Rothman, S. Zaleski, "Lattice-gas models of phase separation: interfaces, phase transitions, and multiphase flow", *Rev. Mod. Phys.* **66** (1994), no. 4, p. 1417.

[12] Y.-L. He, Q. Liu, Q. Li, W.-Q. Tao, "Lattice Boltzmann methods for single-phase and solid-liquid phase-change heat transfer in porous media: a review", *Intl J. Heat Mass Transfer* **129** (2019), p. 160-197.

[13] H. P. G. Darcy, *Les Fontaines publiques de la ville de Dijon. Exposition et application des principes à suivre et des formules à employer dans les questions de distribution d'eau, etc.*, Dalmont, Paris, 1856.

[14] G. Matheron, *Eléments pour une théorie des milieux poreux*, Masson, Paris, 1967.

[15] S. Whitaker, "Flow in porous media I: a theoretical derivation of Darcy's law", *Trans. Porous Med.* **1** (1986), no. 1, p. 3-25.

[16] G. Allaire, "Prolongement de la pression et homogénéisation des équations de Stokes dans un milieu poreux connexe", *C. R. Acad. Sci. Paris* **309** (1989), p. 717-722.

[17] G. Allaire, "One-phase Newtonian flow", in *Homogenization and Porous Media*, Springer, New York, 1997, p. 45-76.

[18] J. Koplik, J. R. Banavar, J. F. Willemsen, "Molecular dynamics of Poiseuille flow and moving contact lines", *Phys. Rev. Lett.* **60** (1988), no. 13, p. 1282.

[19] U. Läcis, P. Johansson, T. Fullana, B. Hess, G. Amberg, S. Bagheri, S. Zaleski, "Steady moving contact line of water over a no-slip substrate: Challenges in benchmarking phase-field and volume-of-fluid methods against molecular dynamics simulations", *Eur. Phys. J. Spec. Top.* **229** (2020), no. 10, p. 1897-1921.

[20] L. Bocquet, J.-L. Barrat, "Flow boundary conditions from nano-to micro-scales", *Soft Matt.* **3** (2007), no. 6, p. 685-693.

[21] C. Cottin-Bizonne, S. Jurine, J. Baudry, J. Crassous, F. Restagno, E. Charlaix, "Nanorheology: an investigation of the boundary condition at hydrophobic and hydrophilic interfaces", *Eur. Phys. J. E* **9** (2002), no. 1, p. 47-53.

[22] L. Bocquet, E. Charlaix, "Nanofluidics, from bulk to interfaces", *Chem. Soc. Rev.* **39** (2010), no. 3, p. 1073-1095.

[23] P. Simonnin, B. Noetinger, C. Nieto-Draghi, V. Marry, B. Rotenberg, "Diffusion under confinement: hydrodynamic finite-size effects in simulation", *J. Chem. Theory Comput.* **13** (2017), no. 6, p. 2881-2889.

[24] R. M. Shroll, D. E. Smith, "Molecular dynamics simulations in the grand canonical ensemble: application to clay mineral swelling", *J. Chem. Phys.* **111** (1999), no. 19, p. 9025-9033.

[25] A. Botan, B. Rotenberg, V. Marry, P. Turq, B. Noetinger, "Carbon dioxide in montmorillonite clay hydrates: thermodynamics, structure, and transport from molecular simulation", *J. Phys. Chem. C* **114** (2010), no. 35, p. 14962-14969.

[26] A. C. Rocha, M. A. Murad, C. Moyne, S. P. Oliveira, T. D. Le, "A new methodology for computing ionic profiles and disjoining pressure in swelling porous media", *Comput. Geosci.* **20** (2016), no. 5, p. 975-996.

[27] T. D. Le, C. Moyne, M. A. Murad, I. Panfilova, "A three-scale poromechanical model for swelling porous media incorporating solvation forces: application to enhanced coalbed methane recovery", *Mech. Mater.* **131** (2019), p. 47-60.

[28] G. Galliéro, J. Colombani, P. A. Bopp, B. Duguay, J.-P. Caltagirone, F. Montel, "Thermal diffusion in micropores by molecular dynamics computer simulations", *Phys. A* **361** (2006), no. 2, p. 494-510.

[29] D. Ameur, G. Galliéro, "Slippage of binary fluid mixtures in a nanopore", *Microfluid. Nanofluid.* **15** (2013), no. 2, p. 183-189.

[30] P. Simonnin, V. Marry, B. Noetinger, C. Nieto-Draghi, B. Rotenberg, "Mineral-and ion-specific effects at clay-water interfaces: structure, diffusion, and hydrodynamics", *J. Phys. Chem. C* **122** (2018), no. 32, p. 18484-18492.
[31] I.-C. Yeh, G. Hummer, "System-size dependence of diffusion coefficients and viscosities from molecular dynamics simulations with periodic boundary conditions", *J. Phys. Chem. B* **108** (2004), no. 40, p. 15873-15879.

[32] O. Plümper, A. Botan, C. Los, Y. Liu, A. Malthe-Sorenssen, B. Jamtveit, "Fluid-driven metamorphism of the continental crust governed by nanoscale fluid flow", *Nat. Geosci.* **10** (2017), no. 9, p. 685-690.

[33] B. Rotenberg, V. Marry, J.-F. Dufrêche, N. Malikova, E. Giffaut, P. Turq, "Modelling water and ion diffusion in clays: a multiscale approach", *C. R. Chim.* **10** (2007), no. 10–11, p. 1108-1116.

[34] R. Shukla, P. Ranjith, A. Haque, X. Choi, "A review of studies on CO₂ sequestration and caprock integrity", *Fuel* **89** (2010), no. 10, p. 2651-2664.

[35] N. Sobecki, C. Nieto-Draghi, A. Di Lella, D. Y. Ding, "Phase behavior of hydrocarbons in nano-pores", *Fluid Phase Equilib.* **497** (2019), p. 104-121.

[36] A. Siria, P. Poncharal, A.-L. Biance, R. Fulcrand, X. Blase, S. T. Purcell, L. Bocquet, "Giant osmotic energy conversion measured in a single transmembrane boron nitride nanotube", *Nature* **494** (2013), no. 7438, p. 455-458.

[37] G. Matheron, "Principles of geostatistics", *Econ. Geol.* **58** (1963), no. 8, p. 1246-1266.

[38] G. De Marsily, F. Delay, J. Gonçalvès, P. Renard, V. Teles, S. Violette, "Dealing with spatial heterogeneity", *Hydrogeol. J.* **13** (2005), no. 1, p. 161-183.

[39] J.-P. Chiles, P. Delfiner, *Geostatistics: Modeling Spatial Uncertainty*, Vol. 497, John Wiley & Sons, New York, 2009.

[40] L. W. Gelhar, *Stochastic Subsurface Hydrology*, Prentice-Hall, New York, 1993.

[41] D. T. Hristopulos, *Random Fields for Spatial Data Modeling A Primer for Scientists and Engineers*, Springer, Netherlands, 2020.

[42] G. Dagan, *Flow and Transport in Porous Formations*, Springer-Verlag GmbH & Co. KG, 1989.

[43] P. Indelman, B. Abramovich, "A higher-order approximation to effective conductivity in media of anisotropic random structure", *Water Resour. Res.* **30** (1994), no. 6, p. 1857-1864.

[44] B. Abramovich, P. Indelman, "Effective permittivity of log-normal isotropic random media", *J. Phys. A* **28** (1995), no. 3, p. 693.

[45] V. V. Jikov, S. M. Kozlov, O. A. Oleinik, *Homogenization of Differential Operators and Integral Functionals*, Springer Science & Business Media, Berlin, Heidelberg, 2012.

[46] S. Armstrong, T. Kuusi, J.-C. Mourrat, *Quantitative Stochastic Homogenization and Large-Scale Regularity*, Vol. 352, Springer, Cham, Switzerland, 2019.

[47] P. King, "The use of renormalization for calculating effective permeability", *Trans. Porous Med.* **4** (1989), no. 1, p. 37-58.

[48] B. Noetinger, "The effective permeability of a heterogeneous porous medium", *Trans. Porous Med.* **15** (1994), p. 99-127.

[49] D. Hristopulos, G. Christakos, "Renormalization group analysis of permeability upscaling", *Stoch Environ. Res. Risk Assess.* **13** (1999), no. 1-2, p. 131-161.

[50] B. Nøetinger, "Computing the effective permeability of log-normal permeability fields using renormalization methods", *C. R. Acad. Sci. - Ser. IIA - Earth Planet. Sci.* **331** (2000), no. 5, p. 353-357.

[51] S. Attinger, "Generalized coarse graining procedures for flow in porous media", *Comput. Geosci.* **7** (2003), no. 4, p. 253-273.

[52] J. Eberhard, S. Attinger, G. Wittum, "Coarse graining for upscaling of flow in heterogeneous porous media", *Multiscale Model. Simul.* **2** (2004), no. 2, p. 269-301.

[53] E. Teodorovich, "Renormalization group method in the problem of the effective conductivity of a randomly heterogeneous porous medium", *J. Expl Theoret. Phys.* **95** (2002), no. 1, p. 67-76.

[54] Y. A. Stepanyants, E. Teodorovich, "Effective hydraulic conductivity of a randomly heterogeneous porous medium", *Water Resour. Res.* **39** (2003), no. 3, p. 360-373.

[55] J. C. Maxwell, *A Treatise on Electricity and Magnetism*, Vol. 1, Clarendon Press, Oxford, 1873.

[56] L. Landau, E. Lifshitz, *Electrodynamics of Continuous Media*, Vol. 8, Pergamon, New York, 1960, 41-43 pages.

[57] Z. Hashin, S. Shtrikman, "A variational approach to the theory of the effective magnetic permeability of multiphase materials", *J. Appl. Phys.* **33** (1962), no. 10, p. 3125-3131.

[58] B. Berkowitz, I. Balberg, "Percolation theory and its application to groundwater hydrology", *Water Resour. Res.* **29** (1993), no. 4, p. 775-794.

[59] A. Hunt, R. Ewing, B. Ghanbarian, *Percolation Theory for Flow in Porous Media*, Vol. 880, Springer, Cham, Switzerland, 2014.

[60] P. King, "The use of field theoretic methods for the study of flow in a heterogeneous porous medium", *J. Phys. A* **20** (1987), no. 12, p. 3935.

[61] B. Noetinger, Y. Gautier, "Use of the Fourier-Laplace transform and of diagrammatical methods to interpret pumping tests in heterogeneous reservoirs", *Adv. Water Resour.* **21** (1998), no. 7, p. 581-590.

[62] M. Mézard, G. Parisi, M. Virasoro, *Spin Glass Theory and Beyond: An Introduction to the Replica Method and its Applications*, Vol. 9, World Scientific Publishing Company, Singapore, 1987.

[63] I. Colecchio, A. Boschan, A. D. Otero, B. Noetinger, "On the multiscale characterization of effective hydraulic
conductivity in random heterogeneous media: a historical survey and some new perspectives", *Adv. Water Resour.* (2020), article no. 103594.

[64] P. Renard, G. De Marsily, "Calculating equivalent permeability: a review", *Adv. Water Resour.* **20** (1997), no. 5, p. 253-278.

[65] G. Dagan, "Higher-order correction of effective permeability of heterogeneous isotropic formations of lognormal conductivity distribution", *Trans. Porous Med.* **12** (1993), no. 3, p. 279-290.

[66] A. De Wit, "Correlation structure dependence of the effective permeability of heterogeneous porous media", *Phys. Fluids* **7** (1995), no. 11, p. 2553-2562.

[67] S. P. Neuman, S. Orr, "Prediction of steady state flow in nonuniform geologic media by conditional moments: exact nonlocal formalism, effective conductivities, and weak approximation", *Water Resour. Res.* **29** (1993), no. 2, p. 341-364.

[68] D. J. Amit, V. Martin-Mayor, *Field Theory, the Renormalization Group, and Critical Phenomena: Graphs to Computers*, 3rd ed., World Scientific Publishing Company, Singapore, 2005.

[69] E. Charlaix, F. Guyon, S. Roux, "Permeability of a random array of fractures of widely varying apertures", *Trans. Porous Med.* **2** (1987), no. 1, p. 31-43.

[70] D. Stauffer, A. Aharony, *Introduction to Percolation Theory*, revised 2nd ed., Taylor & Francis, London, Philadelphia, 2014.

[71] M. Adda-Bedia, Y. Pomeau, "Crack instabilities of a heated glass strip", *Phys. Rev. E* **52** (1995), no. 4, p. 4105.

[72] G. I. Barenblatt, I. P. Zheltov, I. Kochina, "Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata]", *J. Appl. Math. Mech.* **24** (1960), no. 5, p. 1286-1303.

[73] B. Bourbiaux, "Fractured reservoir simulation: a challenging and rewarding issue", *Oil Gas Sci. Technol.-Revue de l'Institut Français du Pétrole* **65** (2010), no. 2, p. 227-238.

[74] O. Bour, P. Davy, "Connectivity of random fault networks following a power law fault length distribution", *Water Resour. Res.* **33** (1997), no. 7, p. 1567-1583.

[75] J.-R. De Dreuzy, P. Davy, O. Bour, "Hydraulic properties of two-dimensional random fracture networks following power law distributions of length and aperture", *Water Resour. Res.* **38** (2002), no. 12, p. 12-1.

[76] J. Maillot, P. Davy, R. Le Goc, C. Darcel, J.-R. De Dreuzy, "Connectivity, permeability, and channeling in randomly distributed and kinematically defined discrete fracture network models", *Water Resour. Res.* **52** (2016), no. 11, p. 8526-8545.

[77] P. M. Adler, J.-F. Thovert, *Fractures and Fracture Networks*, Vol. 15, Springer Science & Business Media, Dordrecht, 1999.

[78] M. Sahimi, *Flow and Transport in Porous Media and Fractured Rock: From Classical Methods to Modern Approaches*, John Wiley & Sons, Germany, 2011.

[79] A. G. Hunt, M. Sahimi, "Flow, transport, and reaction in porous media: percolation scaling, critical-path analysis, and effective medium approximation", *Rev. Geophys.* **55** (2017), no. 4, p. 993-1078.

[80] J. A. Acuna, Y. C. Yortsos, "Application of fractal geometry to the study of networks of fractures and their pressure transient", *Water Resour. Res.* **31** (1995), no. 3, p. 527-540.

[81] B. Nøetinger, N. Jarrige, "A quasi steady state method for solving transient Darcy flow in complex 3D fractured networks", *J. Comput. Phys.* **231** (2012), no. 1, p. 23-38.

[82] B. Nøetinger, "A quasi steady state method for solving transient Darcy flow in complex 3D fractured networks accounting for matrix to fracture flow", *J. Comput. Phys.* **283** (2015), p. 205-223.

[83] B. Mohar, Y. Alavi, G. Chartrand, O. Oellermann, "The Laplacian spectrum of graphs", *Graph Theory, Combin. Appl.* **2** (1991), no. 871-898, p. 12.

[84] B. Mohar, "Some applications of Laplace eigenvalues of graphs", in *Graph Symmetry*, Springer, Dordrecht, 1997, p. 225-275.

[85] M. Bauer, O. Golinelli, "Core percolation in random graphs: a critical phenomena analysis", *Eur. Phys. J. B* **24** (2001), no. 3, p. 339-352.

[86] M. Bauer, O. Golinelli, "Random incidence matrices: moments of the spectral density", *J. Stat. Phys.* **103** (2001), no. 1-2, p. 301-337.

[87] D. Bauer, L. Talon, A. Ehrlacher, "Computation of the equivalent macroscopic permeability tensor of discrete networks with heterogeneous segment length", *J. Hydraul. Eng.* **134** (2008), no. 6, p. 784-793.

[88] C. Bordenave, M. Lelarge, "Resolvent of large random graphs", *Random Struct. Algorithm.* **37** (2010), no. 3, p. 332-352.

[89] G. Semerjian, L. F. Cugliandolo, "Sparse random matrices: the eigenvalue spectrum revisited", *J. Phys. A* **35** (2002), no. 23, p. 4837.

[90] B. Karrer, M. E. Newman, L. Zdeborová, "Percolation on sparse networks", *Phys. Rev. Lett.* **113** (2014), no. 20, article no. 208702.

[91] S. Karra, D. O'Malley, J. Hyman, H. S. Viswanathan, G. Srinivasan, "Modeling flow and transport in fracture networks using graphs", *Phys. Rev. E* **97** (2018), no. 3, article no. 033304.
[92] J. D. Hyman, A. Hagberg, G. Srinivasan, J. Mohd-Yusof, H. Viswanathan, "Predictions of first passage times in sparse discrete fracture networks using graph-based reductions", Phys. Rev. E 96 (2017), no. 1, article no. 013304.

[93] J. D. Hyman, M. Dentz, A. Hagberg, P. K. Kang, "Emergence of stable laws for first passage times in three-dimensional random fracture networks", Phys. Rev. Lett. 123 (2019), no. 24, article no. 248501.

[94] P. Landereau, B. Noetinger, M. Quintard, "Quasi-steady two-equation models for diffusive transport in fractured porous media: large-scale properties for densely fractured systems", Adv. Water Resour. 24 (2001), no. 8, p. 863-876.

[95] R. Haggerty, S. M. Gorelick, "Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity", Water Resour. Res. 31 (1995), no. 10, p. 2383-2400.

[96] T. Babey, J.-R. De Dreuzy, C. Casenave, "Multi-rate mass transfer (MRMT) models for general diffusive porosity structures", Adv. Water Resour. 76 (2015), p. 146-156.

[97] B. Noetinger, T. Estebenet, "Up-scaling of double porosity fractured media using continuous-time random walks methods", Trans. Porous Med. 39 (2000), no. 3, p. 315-337.

[98] B. Noetinger, T. Estebenet, P. Landereau, "A direct determination of the transient exchange term of fractured media using a continuous time random walk method", Trans. Porous Med. 44 (2001), no. 3, p. 539-557.

[99] B. Noetinger, D. Roubinet, A. Russian, T. Le Borgne, F. Delay, M. Dentz, J.-R. De Dreuzy, P. Gouze, "Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale", Trans. Porous Med. 115 (2016), no. 2, p. 345-385.

[100] M. L. Mehta, Random Matrices, Elsevier, Netherlands, 2004.

[101] T. Rogers, I. P. Castillo, R. Kühn, K. Takeda, "Cavity approach to the spectral density of sparse symmetric random matrices", Phys. Rev. E 78 (2008), no. 3, article no. 031116.

[102] M. Biskup et al., "Recent progress on the random conductance model", Probab. Surv. 8 (2011), p. 294-373.

[103] M. Potters, J.-P. Bouchaud, A First Course in Random Matrix Theory, Cambridge University Press, 2019, in press.

[104] M. Valera, Z. Guo, P. Kelly, S. Matz, V. A. Cantu, A. G. Percus, J. D. Hyman, G. Srinivasan, H. S. Viswanathan, "Machine learning for graph-based representations of three-dimensional discrete fracture networks", Comput. Geosci. 22 (2018), no. 3, p. 695-710.

[105] D. O'Malley, S. Karra, J. Hyman, H. S. Viswanathan, G. Srinivasan, "Efficient Monte Carlo with graph-based subsurface flow and transport models", Water Resour. Res. 54 (2018), no. 5, p. 3758-3766.

[106] R. Romeu, B. Noetinger, "Calculation of internodal transmissivities in finite difference models of flow in heterogeneous porous media", Water Resour. Res. 31 (1995), no. 4, p. 943-959.

[107] G. Biroli, J.-P. Bouchaud, M. Potters, "Extreme value problems in random matrix theory and other disordered systems", J. Statist. Mech.: Theory Exp. 2007 (2007), no. 07, article no. 07019.

[108] V. A. Marchenko, L. A. Pastur, "Distribution of eigenvalues for some sets of random matrices", Mat. Sborn. 114 (1967), no. 4, p. 507-536.

[109] C. Louart, Z. Liao, R. Couillet et al., "A random matrix approach to neural networks", Ann. Appl. Probab. 28 (2018), no. 2, p. 1190-1248.

[110] L. Dall'Amico, R. Couillet, N. Tremblay, "Classification spectrale par la laplacienne déformée dans des graphes réalistes", in XXVIᵉ colloque GRETSI (GRETSI 2019), Aug 2019, lille, France, 2019, hal-02153901.

[111] D. L. Koch, J. F. Brady, "Dispersion in fixed beds", J. Fluid Mech. 154 (1985), p. 399-427.

[112] D. L. Koch, J. F. Brady, "Anomalous diffusion in heterogeneous porous media", Phys. Fluids 31 (1988), no. 5, p. 965-973.

[113] G. I. Taylor, "Conditions under which dispersion of a solute in a stream of solvent can be used to measure molecular diffusion", Proc. R. Soc. Lond. A. Math. Phys. Sci. 225 (1954), no. 1163, p. 473-477.

[114] P. Saffman, "Dispersion due to molecular diffusion and macroscopic mixing in flow through a network of capillaries", J. Fluid Mech. 7 (1960), no. 2, p. 194-208.

[115] C. Baudet, E. Guyon, Y. Pomeau, "Dispersion dans un écoulement de Stokes", J. Phys. Lett. 46 (1985), no. 21, p. 991-998.

[116] E. Flekkøy, U. Oxaal, J. Feder, T. Jøssang, "Hydrodynamic dispersion at stagnation points: simulations and experiments", Phys. Rev. E 52 (1995), no. 5, p. 4952.

[117] F. Gjetvaj, A. Russian, P. Gouze, M. Dentz, "Dual control of flow field heterogeneity and immobile porosity on non-Fickian transport in Berea sandstone", Water Resour. Res. 51 (2015), no. 10, p. 8273-8293.

[118] G. Matheron, G. De Marsily, "Is transport in porous media always diffusive? A counterexample", Water Resour. Res. 16 (1980), no. 5, p. 901-917.

[119] J.-P. Bouchaud, A. Georges, "Anomalous diffusion in disordered media: statistical mechanisms, models and physical applications", Phys. Rep. 195 (1990), no. 4-5, p. 127-293.

[120] J.-P. Bouchaud, A. Georges, J. Koplik, A. Provata, S. Redner, "Superdiffusion in random velocity fields", Phys. Rev. Lett. 64 (1990), no. 21, p. 2503.

[121] U. Jaekel, H. Vereecken, "Renormalization group analysis of macrodispersion in a directed random flow", Water Resour. Res. 33 (1997), no. 10, p. 2287-2299.
[122] M. D. Hürlimann, L. M. Schwartz, P. N. Sen, "Probability of return to the origin at short times: a probe of microstructure in porous media", *Phys. Rev. B* **51** (1995), no. 21, p. 14936.

[123] N. Krepysheva, L. Di Pietro, M.-C. Néel, "Space-fractional advection–diffusion and reflective boundary condition", *Phys. Rev. E* **73** (2006), no. 2, article no. 021104.

[124] A. Zoia, M.-C. Néel, A. Cortis, "Continuous-time random-walk model of transport in variably saturated heterogeneous porous media", *Phys. Rev. E* **81** (2010), no. 3, article no. 031104.

[125] V. Guillon, M. Fleury, D. Bauer, M.-C. Neel, "Superdiffusion in homogeneous unsaturated porous media using NMR propagators", *Phys. Rev. E* **87** (2013), no. 4, article no. 043007.

[126] M.-C. Néel, D. Bauer, M. Fleury, "Model to interpret pulsed-field-gradient NMR data including memory and superdispersion effects", *Phys. Rev. E* **89** (2014), no. 6, article no. 062121.

[127] M. Dentz, T. Le Borgne, A. Englert, B. Bijeljic, "Mixing, spreading and reaction in heterogeneous media: a brief review", *J. Contam. Hydrol.* **120** (2011), p. 1-17.

[128] T. Le Borgne, M. Dentz, D. Bolster, J. Carrera, J.-R. de Dreuzy, P. Davy, "Non-Fickian mixing: temporal evolution of the scalar dissipation rate in heterogeneous porous media", *Adv. Water Resour.* **33** (2010), no. 12, p. 1468-1475.

[129] T. Le Borgne, M. Dentz, J. Carrera, "Lagrangian statistical model for transport in highly heterogeneous velocity fields", *Phys. Rev. Lett.* **101** (2008), no. 9, article no. 090601.

[130] M. Dentz, P. K. Kang, A. Comolli, T. Le Borgne, D. R. Lester, "Continuous time random walks for the evolution of Lagrangian velocities", *Phys. Rev. Fluids* **1** (2016), no. 7, article no. 074004.

[131] T. Le Borgne, P. D. Huck, M. Dentz, E. Villermaux, "Scalar gradients in stirred mixtures and the deconstruction of random fields", *J. Fluid Mech.* **812** (2017), p. 578-610.

[132] P. G. Saffman, G. I. Taylor, "The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid", *Proc. R. Soc. Lond. A* **245** (1958), no. 1242, p. 312-329.

[133] G. M. Homsy, "Viscous fingering in porous media", *Annu. Rev. Fluid Mech.* **19** (1987), no. 1, p. 271-311.

[134] C. Tang, "Diffusion-limited aggregation and the Saffman–Taylor problem", *Phys. Rev. A* **31** (1985), no. 3, p. 1977.

[135] B. I. Shraiman, "Velocity selection and the Saffman–Taylor problem", *Phys. Rev. Lett.* **56** (1986), no. 19, p. 2028.

[136] J. Langer, "Dendrites, viscous fingers, and the theory of pattern formation", *Science* **243** (1989), no. 4895, p. 1150-1156.

[137] M. King, V. A. Dunayevsky, "Why waterflood works: a linearized stability analysis", in *SPE Annual Technical Conference and Exhibition*, Society of Petroleum Engineers, 1989.

[138] A. De Wit, G. Homsy, "Viscous fingering in periodically heterogeneous porous media. I. Formulation and linear instability", *J. Chem. Phys.* **107** (1997), no. 22, p. 9609-9618.

[139] A. De Wit, G. Homsy, "Viscous fingering in periodically heterogeneous porous media. II. Numerical simulations", *J. Chem. Phys.* **107** (1997), no. 22, p. 9619-9628.

[140] V. Artus, B. Nøetinger, L. Ricard, "Dynamics of the water–oil front for two-phase, immiscible flow in heterogeneous porous media. 1–stratified media", *Trans. Porous Med.* **56** (2004), no. 3, p. 283-303.

[141] B. Nøetinger, V. Artus, L. Ricard, "Dynamics of the water–oil front for two-phase, immiscible flow in heterogeneous porous media. 2–isotropic media", *Trans. Porous Med.* **56** (2004), no. 3, p. 305-328.

[142] V. Artus, B. Noetinger, "Up-scaling two-phase flow in heterogeneous reservoirs: current trends", *Oil Gas Sci. Technol.* **59** (2004), no. 2, p. 185-195.

[143] K. T. Tallakstad, H. A. Knudsen, T. Ramstad, G. Løvoll, K. J. Måløy, R. Toussaint, E. G. Flekkøy, "Steady-state two-phase flow in porous media: statistics and transport properties", *Phys. Rev. Lett.* **102** (2009), no. 7, article no. 074502.

[144] R. Toussaint, K. J. Måløy, Y. Méheust, G. Løvoll, M. Jankov, G. Schäfer, J. Schmittbuhl, "Two-phase flow: structure, upscaling, and consequences for macroscopic transport properties", *Vadose Zone J.* **11** (2012), no. 3, article no. 2011.0123.

[145] E. Koval et al., "A method for predicting the performance of unstable miscible displacement in heterogeneous media", *Soc. Petrol. Eng. J.* **3** (1963), no. 02, p. 145-154.

[146] M. Todd, W. Longstaff et al., "The development, testing, and application of a numerical simulator for predicting miscible flood performance", *J. Petrol. Tech.* **24** (1972), no. 07, p. 874-882.

[147] Y. C. Yortsos, "A theoretical analysis of vertical flow equilibrium", *Trans. Porous Med.* **18** (1995), no. 2, p. 107-129.

[148] M. Blunt, M. Christie, "How to predict viscous fingering in three component flow", *Trans. Porous Med.* **12** (1993), no. 3, p. 207-236.

[149] K. Sorbie, H. Zhang, N. Tsibuklis, "Linear viscous fingering: new experimental results, direct simulation and the evaluation of averaged models", *Chem. Eng. Sci.* **50** (1995), no. 4, p. 601-616.

[150] T. Witten Jr, L. M. Sander, "Diffusion-limited aggregation, a kinetic critical phenomenon", *Phys. Rev. Lett.* **47** (1981), no. 19, p. 1400.

[151] D. Wilkinson, J. F. Willemse, "Invasion percolation: a new form of percolation theory", *J. Phys. A* **16** (1983), no. 14, p. 3365.

[152] L. Paterson, "Diffusion-limited aggregation and two-fluid displacements in porous media", *Phys. Rev. Lett.* **52** (1984), no. 18, p. 1621.
[153] J. G. Masek, D. L. Turcotte, "A diffusion-limited aggregation model for the evolution of drainage networks", *Earth Planet. Sci. Lett.* **119** (1993), no. 3, p. 379-386.

[154] S. Saha, S. Atis, D. Salin, L. Talon, "Phase diagram of sustained wave fronts opposing the flow in disordered porous media", *Europhys. Lett.* **101** (2013), no. 3, p. 38003.

[155] S. Atis, A. K. Dubey, D. Salin, L. Talon, P. Le Doussal, K. J. Wiese, "Experimental evidence for three universality classes for reaction fronts in disordered flows", *Phys. Rev. Lett.* **114** (2015), no. 23, article no. 234502.

[156] M. Kardar, G. Parisi, Y.-C. Zhang, "Dynamic scaling of growing interfaces", *Phys. Rev. Lett.* **56** (1986), no. 9, p. 889.

[157] B. Noetinger, G. Zargar, "Multiscale description and upscaling of fluid flow in subsurface reservoirs", *Oil Gas Sci. Technol.* **59** (2004), no. 2, p. 119-139.

[158] C. E. Cohen, D. Ding, M. Quintard, B. Bazin, "From pore scale to wellbore scale: impact of geometry on wormhole growth in carbonate acidization", *Chem. Eng. Sci.* **63** (2008), no. 12, p. 3088-3099.

[159] A. De Wit, "Chemo-hydrodynamic patterns in porous media", *Phil. Trans. R. Soc. A* **374** (2016), no. 2078, article no. 20150419.

[160] L. de Arcangelis, S. Redner, H. Herrmann, "A random fuse model for breaking processes", *J. Phys. Lett.* **46** (1985), no. 13, p. 585-590.

[161] L. de Arcangelis, H. Herrmann, "Scaling and multiscaling laws in random fuse networks", *Phys. Rev. B* **39** (1989), no. 4, p. 2678.

[162] P. L. Krapivsky, S. Redner, F. Leyvraz, "Connectivity of growing random networks", *Phys. Rev. Lett.* **85** (2000), no. 21, p. 4629.

[163] P. S. Dodds, D. H. Rothman, "Scaling, universality, and geomorphology", *Annu. Rev. Earth Planet. Sci.* **28** (2000), no. 1, p. 571-610.

[164] D. L. Turcotte, "Self-organized complexity in geomorphology: observations and models", *Geomorphology* **91** (2007), no. 3-4, p. 302-310.

[165] M. Keiler, J. Knight, S. Harrison, "Climate change and geomorphological hazards in the eastern European Alps", *Phil. Trans. R. Soc. A* **368** (2010), no. 1919, p. 2461-2479.

[166] A. Tarantola, *Inverse Problem Theory and Methods for Model Parameter Estimation*, Vol. 89, SIAM, Philadelphia, 2005.

[167] A. M. Lavenue, B. S. Ramarao, G. De Marsily, M. G. Marietta, "Pilot point methodology for automated calibration of an ensemble of conditionally simulated transmissivity fields: 2. Application", *Water Resour. Res.* **31** (1995), no. 3, p. 495-516.

[168] A. Abellan, B. Noetinger, "Optimizing subsurface field data acquisition using information theory", *Math. Geosci.* **42** (2010), no. 6, p. 603-630.

[169] L. Zdeborová, F. Krzakala, "Statistical physics of inference: thresholds and algorithms", *Adv. Phys.* **65** (2016), no. 5, p. 453-552.

[170] L. Zdeborová, F. Krząkała, "Phase transitions in the coloring of random graphs", *Phys. Rev. E* **76** (2007), no. 3, article no. 031131.