

2021/2023 UNC Greensboro PDE Conference,

Electronic Journal of Differential Equations, Conference 26 (2022/2025), pp. 179–200.

ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu

DOI: 10.58997/ejde.conf.26.g2

p-LAPLACIAN IN PHENOMENOLOGICAL MODELING OF

FLOW IN POROUS MEDIA AND CFD SIMULATIONS

PETR GIRG, LUKÁŠ KOTRLA, ANEŽKA ŠVANDOVÁ

Abstract. The aim of this article is to discuss several aspects of connections

between the p-Laplacian and mathematical models in hydrology. At first we
present models of groundwater flow in phreatic aquifers and models of irriga-

tion and drainage that lead to quasilinear parabolic equations involving the

p-Laplacian. Next, we survey conditions of validity of Strong Maximum Prin-
ciple and Strong Comparison Principle for this type of problems. Finally, we

employ computer fluid dynamics simulations to realistic scenario of fracture

networks to estimate values of the parameters of constitutive laws governing
groundwater flow in the context of fractured hard-rock aquifers.

1. Introduction

The aim of this article is to discuss several aspects of connections between the
p-Laplacian and mathematical models of groundwater flow with applications to
irrigation, drainage, and fresh water supply in small rural areas. The operator
p-Laplacian, a non-linear generalization of the Laplace operator,

u 7→ div
(
|∇u|p−2∇u

)
, p > 1, p ̸= 2,

gained substantial attention of mathematicians working in nonlinear functional
analysis since the 1970s. This can be attributed mostly to the fact that the p-
Laplacian exhibits many features of more general nonlinear operators, while being
simple and elegant, and hence allowing the ideas of the proofs to be kept clear and
accessible. However, origins of the p-Laplacian can be traced back to the 1870s,
see, e.g., [54, 55, 61], when the radially symmetric version of this operator (written
in a different way than it is customary today) was used in the theoretical study of
groundwater flow towards a well in coarse grained porous media (such as gravels).
To the best of our knowledge, the p-Laplacian in its full PDE form was first intro-
duced in the paper [34] in the context of mathematical models of the flow of natural
gas in the porous rock as early as in 1945. In particular, the following equation
(written in modern notation)

u−κ(x, t)
∂u

∂t
(x, t) = c∆p u(x, t) ,

2020 Mathematics Subject Classification. 76S05, 35Q35, 35K92.
Key words and phrases. p-Laplacian; porous medium; filtration; Darcy’s law;
pressure-to-velocity power law.
©2025 This work is licensed under a CC BY 4.0 license.

Published May 13, 2025.
179


180 P. GIRG, L. KOTRLA, A. ŠVANDOVÁ EJDE-2022/25/CONF/26

where c > 0, κ = 1/2, p = 3/2, x ∈ R3, and t ∈ (0,+∞), was proposed as a model
of isothermic process of turbulent filtration of natural gas in a porous medium in
[34]. This equation was soon after suggested in a more general setting as a model of
polytropic process of turbulent filtration of natural gas with κ = γ/(γ + 1), where
γ is the polytropic index of the process and p ∈ [3/2, 2], see [35, p. 504]. In the
early papers on the subject, authors (without the benefit of techniques of modern
functional analysis) limited themselves to study particular cases of initial-boundary
value problems with the p-Laplacian which were important for engineers of that
days, with breakthrough results in [4, 5]. More information about pioneering works
related to applications of the p-Laplacian to flow in porous media can be found in
[11].

In this article, we cover two topics connecting the p-Laplacian and mathematical
models in hydrology. The first topic is related to mathematical models in irrigation
and drainage and concerns validity of Strong Maximum Principle (SMP) and Strong
Comparison Principle (SCP) for initial-value problem

∂

∂t
b(u(x, t))− c∆pu(x, t) = f(x, t) ≥ 0 for (x, t) ∈ Ω× (0, T );

u(x, 0) = u0(x) ≥ 0 for x ∈ Ω;

u(x, t) = 0 for (x, t) ∈ ∂Ω× (0, T ) ,

u(x, t) ≥ 0 for (x, t) ∈ Ω× [0, T ] ,

(1.1)

where Ω ⊂ RN , N ∈ N, is a bounded domain, c ≡ const. > 0, b : R+ → R+ :=
[0,+∞) is a continuous function, b(0) = 0, and b ∈ C1(0,+∞) with b′ > 0 in
(0,+∞), f ∈ C(Ω × (0, T )) and u0 ∈ C(Ω) are nonnegative functions. For a de-
tailed discussion about mathematical models used in irrigation and drainage based
on (1.1), see Section 2. Explicit solutions of problem (1.1) can be obtained only
in rare special cases. Thus we rely on qualitative analysis and numerical methods
in dealing with (1.1). Maximum and comparison principles play important role
in this process. In this paper, we discuss conditions under which these principles
for (1.1) hold and provide some realistic examples when they do not hold. This is
not interesting only from theoretical point of view, but it has implications for the
choices of appropriate numerical methods to find numerical solutions of (1.1).

The second topic delves into the utilization of computer fluid dynamics (CFD)
numerical simulations to estimate values of the parameters of constitutive laws
governing groundwater flow in the context of fractured hard rock. This research
builds upon existing research, see, e.g., [15, 47, 57, 60, 63], and is motivated by
the fact that increasing demand for fresh water has driven interest in hard-rock
aquifers [24, 50], despite their limited well yields due to water flow occurring only
in cracks and fractures. With hard-rocks covering over 20% of the Earth’s landmass
[3, 31, 24], these aquifers hold significant freshwater reserves, particularly in semi-
arid regions like sub-Saharan Africa [38, 59], Australia [22] and India [44, 53]. They
provide important water sources for rural populations in these areas, particularly
in Australia and India, which store 40% and 50% of their groundwater in such
aquifers, respectively. Understanding flow dynamics in these aquifers is essential
for improving rural living standards through better access to fresh water.

This article is organized as follows. In Section 2, we present several mathemat-
ical models of groundwater flow in phreatic aquifers and related models used in


EJDE-2022/2025/CONF/26 FLOW IN POROUS MEDIA AND CFD SIMULATIONS 181

irrigation and drainage. In Section 3, we survey some of our recent results concern-
ing SMP and SCP. In Section 4, we present some of our recent results concerning
estimations of parameters in constitutive laws governing groundwater flow through
fractured hard rocks. Section 5 is devoted to contribution of P. Drábek to topics
related to the p-Laplacian. This paper adheres to the SI system for all physical
quantities (m for length, s for time, kg for mass, etc.).

2. Mathematical models

In mathematical models of groundwater flow, averaged velocity, total head, and
piezometric head are crucial concepts. Actual velocity of the groundwater highly
oscillates in the channels in the porous medium and thus it is difficult to mea-
sure and predict. Averaged velocity v⃗av captures the idea of bulk motion of the
groundwater within a sufficiently large control volume in the porous medium. It is
defined by means of specific discharge q⃗ by formula v⃗av = q⃗/ϕeff , where the specific
discharge (vector) q⃗ takes the direction of the flow and its magnitude is defined
as the volume of water flowing per unit time through a unit cross-sectional area
normal to the direction of flow, and ϕeff is effective porosity of the medium see,
e.g., [7, p. 121] for detailed explanation. On this macroscopic level, the process
of transformation and dissipation of energy can be described in the following way.
The total mechanical energy per volume ET in a control volume of water is the sum
of gravitational potential energy zϱg, pressure energy P , and kinetic energy 1

2ϱv
2
av

(all three per volume), where vav stands for the magnitude of averaged velocity of
the flow in the control volume, ϱ is the water density, P pressure, z elevation of
the control volume from the datum, g gravitational acceleration, see, e.g., [49]. For
incompressible liquid such as water, one can equivalently consider another quantity

hT =
ET

ϱg
= z +

P

ϱg
+

1

2g
v2av

which can be directly measured in practice, e.g., by using observation wells or by so
called piezometers, see, e.g., [7, p. 63] or [49, pp. 129–132]. Groundwater is loosing
its total energy (or equivalently total head hT ) while flowing due to viscous forces
and friction with porous medium. Thus, its total energy decreases in the direction
of the flow. In typical real-world situations, the term corresponding to kinetic
energy is negligible and can be dropped, see, e.g., [25, p. 5] or [49, pp. 40–43]. In
this way, we obtain piezometric head

h = z +
P

ϱg

which is the state variable in the mathematical models of the water flow in the
underground. On the other hand, the specific discharge is the flux quantity. The
constitutive law relating this two quantities, quantitatively describes the rate of
dissipation of the energy along the flow path.

In this paper, we will limit ourselves to a mathematical model of unconfined
aquifer bounded below by a flat impermeable layer at z = 0. This constitutes a free
boundary problem in its full generality, since the upper boundary is the unknown
surface of the groundwater. In groundwater modelling, this difficulty is simplified by
assuming that the vertical flux in the aquifer is negligible, leading to piezometric
head h being constant in the z-direction. This assumption is known as Dupuit-
Forchheimer assumption, see, e.g., [7, 8, 9] and it is based on observations performed


182 P. GIRG, L. KOTRLA, A. ŠVANDOVÁ EJDE-2022/25/CONF/26

on aquifers. Then the height ĥ of the free surface of the groundwater (called water

table for short) above the point (x, y, 0) at time t is ĥ(x, y, t) = h(x, y, 0, t). The
balance equation can be written as follows

ϕeff
∂ĥ

∂t
(x, y, t)− div

(
ĥ (x, y, t) q⃗ (x, y, t)

)
= ĝ(x, y, t) , (2.1)

see, e.g., [8, Eq. (5.4.43)], [9, Eq. 2.6]. Here, ĝ(x, y, t) represents external sinks
and sources (evaporation, rainfall etc.) of volume of water per area and time unit.
Recall that ϕeff is effective porosity of the porous medium and q⃗(x, y, t) is specific
discharge.

To eliminate the flux variable from the balance equation (2.1), suitable consti-
tutive law is needed. It is obtained empirically from experimental data for given
porous medium and fluid. In practice, these experiments are performed on a sample
of porous medium subjected to one dimensional flow for several values of magnitudes
of flux q, which are kept constant during each measurement. Linear Darcy’s law is
the most widely used in practice due to its simplicity and still reasonable accuracy.
It relates groundwater flux to the piezometric head loss per length according to the
following formula

q = c
△h
△L

, (2.2)

where c > 0 is a constant to be determined from measured data, △h is the difference
of the piezometric heads measured at two distinct locations distance △L apart.
This formula was established experimentally for filtration of water through sand by
Henry Darcy [19] in 1856. Later, it was found that it has limited range of its validity
in coarse grained media (such as gravels), see, e.g., [20, 26, 30, 41, 54, 55] as well as in
media with very low permeability (such as clays, certain soils, and sandstones), see,
e.g., [17, 28, 64]. For thorough surveys and discussions of this and other constitutive
laws and various criteria of their validity, see, e.g., [1, 7, 48, 49, 56]. It follows from
discussions in these papers that the power-type law

q = c
(△h
△L

)m

(2.3)

with constants c,m > 0 to be determined from measured data, is simple but flexible
enough to fit with most experimental data obtained for various porous media. For
m = 1, the power-type law coincides with Darcy’s law. The case m > 1 correponds
to natural media such as clays, certain soils and sandstones, while the case 1/2 <
m < 1 corresponds to coarse grained materials such as gravels, see the literature
listed above. Let us note that the law (2.3) is not the only type of nonlinear laws
used in practice. For thorough surveys, see, e.g., [1, 11, 7, 48, 49].

For the homogeneous and isotropic porous medium, the two- or three-dimensional
constitutive law in differential form can be inferred from the one-dimensional one,
by taking into account that the flux takes the opposite direction of the gradient of
the piezometric head and no flow occurs if the gradient of the piezometric head is
zero. In this way, we obtain

q⃗ =

{
0⃗ for ∇h = 0⃗ ,

−c |∇h|m−1 ∇h for ∇h ̸= 0⃗ ,
(2.4)

where ∇h stands for the spatial gradient of the piezometric head, c,m > 0 are
constants as in (2.3).


EJDE-2022/2025/CONF/26 FLOW IN POROUS MEDIA AND CFD SIMULATIONS 183

By substituting two-dimensional differential version of Darcy’s law, i.e., (2.4)
with m = 1 for q⃗ in (2.1), we obtain classical porous medium equation

ϕeff
∂ĥ

∂t
− cdiv(ĥ∇ĥ) = ĝ(x, y, t) .

By the same process, we obtain the nonlinear equation

ϕeff
∂ĥ

∂t
− cdiv(ĥ|∇ĥ|m−1∇ĥ) = ĝ(x, y, t) , (2.5)

by using (2.4) with m > 0, m ̸= 1.
Now we turn our attention to mathematical models from irrigation and drainage.

We were motivated by models proposed in [46], but we take into account nonlinear
effects and use the power-type law (2.3) instead of Darcy’s law (2.2).

At first, let us consider local aquifer under a field, represented by a bounded
domain Ω ⊂ R2, surrounded by open water body as in Figure 1. The horizon-
tal permeable layer is bounded from below by the horizontal impermeable layer
(bedrock). The drainage channels are fully penetrating, i.e., they reach the surface
of the impermeable bedrock. We take the bedrock surface as vertical datum, i.e., we
assign to it coordinate z = 0. Let H > 0 be the depth of the open water bodies rel-

ative to the datum. Then the height ĥ of the water table above the datum satisfies

balance equation (2.5) with boundary conditions ĥ = H on ∂Ω. Using substitution

p = m+ 1 (to match notation with the p-Laplacian), u = ĥp/(p−1) −Hp/(p−1), we
obtain initial-boundary value problem (1.1) with

b(u) = ϕeff

( p

p− 1

)p−1[(
u+H

p
p−1

) p−1
p −H

]
, (2.6)

f(x, t) =
(

p
p−1

)p−1
ĝ(x, t), and u0(x) = ĥ

p/(p−1)
0 (x)−Hp/(p−1), where x = (x, y) ∈

Ω.
The second model covered by our theoretical methods concerns local phreatic

aquifer under a strip field between two parallel ditches. The aquifer is bounded
from below by an impermeable bedrock and the ditches are fully penetrating to
the bedrock, see Figure 2. In general, the domain would be infinite strip in this
situation. However, if we assume translation invariance of the solution along the
parallel ditches, we can limit ourselves to a bounded interval Ω = (−L/2, L/2),
where L > 0, see Figure 2. This nonlinear model has been proposed in [9] and it is

motivated by [6, 39, 51, 52], where linear Darcy’s law (2.2) is used instead. With ĥ

being piezometric head and the bedrock vertical datum, the function u = ĥp/(p−1)−
Hp/(p−1) satisfies (1.1) with b, f, u0 as above but with x = x ∈ Ω = (−L/2, L/2).

3. Maximum and comparison principles

We will address the question of validity of SMP for a nonnegative weak solution
to

∂

∂t
b(u(x, t))−∆pu(x, t) = f(x, t) for (x, t) ∈ Ω× (0, T );

u(x, 0) = u0(x) for x ∈ Ω;

u(x, t) = 0 for (x, t) ∈ ∂Ω× (0, T ),

which is (1.1) with c = 1 for simplicity. The reader is referred, e.g., to [9] for
the precise definition of weak solution and the appropriate function spaces. Recall


184 P. GIRG, L. KOTRLA, A. ŠVANDOVÁ EJDE-2022/25/CONF/26

(a)

Ω

(b)

Figure 1. Fields (brown) drained by a grid of ditches (blue). (a)
Topographical view. Domain Ω represents a field surrounded by a
free water body. (b) 3D detail of a field surrounded by free water
body accumulated in ditches.

that we assume the following hypotheses: b : R+ → R+ is a continuous function,
b(0) = 0, and b ∈ C1(0,+∞) with b′ > 0 in (0,+∞). For simplicity, we assume that
both, f : Ω × (0, T ) → R and u0 : Ω → R, are continuous and nonnegative. Then
indeed, u ≥ 0 on Ω× (0, T ) by Weak Comparison Principle (WCP) [9, Proposition
1, p. 28].

In the modelling of groundwater flow, the way the water spreads from the source
area (where rain or irrigation occurs) is crucial question. In particular, we are
interested in the following scenario. The initial distribution of water u0 ∈ C (Ω) is
such that u0 > 0 on some compact subset of Ω. Considering diffusion process, one
would expect that the groundwater immediately spreads toward the boundary ∂Ω.
Indeed, the solution to linear diffusion equation satisfies the following SMP.


EJDE-2022/2025/CONF/26 FLOW IN POROUS MEDIA AND CFD SIMULATIONS 185

(a)

(b)

J

H Hh

Rain

Water TableDitch Ditch

Impermeable

BedrockL
Ω=(-L/2, L/2)

Figure 2. Fields (brown) drained by parallel ditches (blue). (a)
Topographical view. (b) 2D detail of a field between two ditches.

Definition 3.1 (Strong maximum principle). Let f ≥ 0 on Ω × (0, T ) and u ∈
C
(
Ω× [0, T )

)
be a corresponding nonnegative weak solution to (1.1). We say that

u satisfies SMP if there exists τ ∈ (0, T ) such that u(x, t) > 0 for all (x, t) ∈ Ω×(0, τ)
and u(x, t) ≡ 0 for all (x, t) ∈ Ω× [τ, T ).

In the nonlinear case p ̸= 2 and/or b(s) ̸≡ s, the situation is not so clear-cut at
all. The validity of SMP depends on p and properties of b(s) as s → 0+. At first,
we recall famous Barenblatt’s self-similar solution of Leibenson’s equation obtained
in [4].

Example 3.2 (Barenblatt’s self-similar solution). Let N ∈ N, k, p > 1 be given
constants. Existence of self-similar radially symmetric solution of Leibenson’s equa-
tion

∂u1/k

∂t
= ∆pu in RN × (0,+∞) ,

was established in [4]. In particular, it was shown using ODE techniques that the
following radially symmetric problem

∂

∂t
w1/k =

1

rN−1

∂

∂r

[
rN−1

∣∣∂w
∂r

∣∣p−2 ∂w

∂r

]
in (0,+∞)× (0,+∞) ,


186 P. GIRG, L. KOTRLA, A. ŠVANDOVÁ EJDE-2022/25/CONF/26

possesses self-similar solution, i.e., solution of the form

w(r, t) = t−
N
β BN,k,p

( r

t1/β

)
,

where β = β(N, k, p) > 0, which satisfies
∫ +∞
0

w(r, t)rN−1dr = 1 for any t > 0.
Then, for a self-similar radially symmetric solution of Leibenson’s equation, it holds
u(x, t) = w(|x|, t) for x ∈ RN and t > 0. Moreover, it easily follows that the initial
trace of u is the Dirac measure in RNconcentrated at 0, i.e., limt→0+ u(·, t) → δ0
in the sense of measures in RN .

Three qualitatively different cases need to be distinguished here.

(i) If k(p− 1) > 1, then

BN,k,p(s) =
[
max

{
0, C − κs

p
p−1

}]γ
,

where γ = γ(N, k, p) > 0, κ = κ(N, k, p) > 0, C = C(N, k, p) > 0.
(ii) If k(p− 1) = 1, then β = p and

BN,k,p(s) = C exp
(
−

(s
p

) p
p−1

)
,

where C = C(N, k, p) > 0.
(iii) If 0 < k(p− 1) < 1, then

BN,k,p(s) =
(
C + κs

p
p−1

)γ

,

where γ = γ(N, k, p) > 0, κ = κ(N, k, p) > 0, C = C(N, k, p) > 0.

It can be easily seen that the self-similar solution w has compact support suppw =

[0, t
1
β
(
C
κ

) p−1
p ] in the case (i), while w(r, t) > 0 for any r ≥ 0, t > 0, in the cases (ii)

and (iii). This makes the case (i) qualitatively distinct from the other two cases.
Now we turn our eyes back to (1.1). Let k, p > 1, k(p− 1) > 1, b(s) = s1/k, and

x0 ∈ Ω be arbitrary but fixed. Let also 0 < σ < 1
2 (κ/C)

β(p−1)/p
dist(x0, ∂Ω)

β
be

fixed. Then the function

u(x, t) ≡ w(|x− x0|, t+ σ)

is a solution to (1.1) with f(x, t) ≡ 0 and u0(x) ≡ w(|x− x0|, σ) for all x ∈ Ω and

0 < t < 1
2 (κ/C)

β(p−1)/p
dist(x0, ∂Ω)

β
. Taking into account that u0 ≥ 0 and it is

positive on a set of positive measure, this means that SMP does not hold for (1.1)
for this choice of parameters and function b. This motivated our research in [10, 12]
to obtain sufficient conditions for the validity of SMP as well as discovery of further
counterexamples for SCP.

Indeed, we obtained the following affirmative result for the SMP, see [10] for
more details and the proof.

Proposition 3.3 ([10, Thm 1.1.]). Let 1 < p < 2, N ≥ 1 and assume that b :
R+ → R+ is as above in (1.1) and satisfies

lim
s→0+

s2−pb′(s)

| log s|p−1
= 0. (3.1)

Finally, assume that u : Ω̄ × [0, T ) → R+ is a continuous, nonnegative, weak
solution to (1.1). Then, for any fixed t0 ∈ (0, T ), the solution u (·, t0) is either
positive everywhere on Ω or else identically zero on Ω.


EJDE-2022/2025/CONF/26 FLOW IN POROUS MEDIA AND CFD SIMULATIONS 187

In particular, if u(ξ, 0) = u0(ξ) > 0 for some ξ ∈ Ω, then there exists τ ∈ (0, T ]
such that u(x, t) > 0 for all (x, t) ∈ Ω × (0, τ), i.e., Strong Maximum Principle
is valid in the (N + 1)-dimensional space-time cylinder Ω × (0, τ). The number
τ ∈ (0, T ) can be estimated from below by

τ = sup {T ′ ∈ (0, T ] : u(ξ, t) > 0 for all t ∈ [0, T ′)} > 0 .

Note that b is given by (2.6) in our models. Hence, b′(s) > k ≡ const. for all s ≥ 0
and the validity of SMP depends on value of p only. The following counterexample is
slightly modified [12, Example 2.3, pp. 368–369], where the comparison of stationary
solution u and evolutionary solution v to (1.1) with b(s) ≡ s is studied in one space
dimension. To obtain a counterexample to SMP, we will compare evolutionary
solution with the trivial one.

Example 3.4 (Counterexample to SMP for p > 2 in 1D). Let p > 2, max{2, p/(p−
2)} < β and 0 < γ ≤ β. Assume also that b′(s) ≥ k ≡ const. > 0. In accordance
with [12], define

u(x, t) ≡ u(x) = 0 and

v(x, t) =

{
0 for x ∈ (−1, 0] ,

txβ(1− xγ) for x ∈ (0, 1) .

Clearly, u(±1, t) = v(±1, t) = 0, v(x, 0) = u(x, 0) = 0, u(x, t) = v(x, t) = 0 on
(−1, 0)× (0, T ), and 0 = u(x, t) < v(x, t) on (0, 1)× (0, T ).

Our goal is to show that there exists t0 > 0 such that

f(x, t)
def
= b′(v)

∂v

∂t
− ∂

∂x

(∣∣∂v
∂x

∣∣p−2 ∂v

∂x

)
> 0 on (0, 1)× (0, t0) .

Indeed,

f(x, t) ≥ k
∂v

∂t
− ∂

∂x

(∣∣∂v
∂x

∣∣p−2 ∂v

∂x

)
= kxβ

[
(1− xγ)− (p− 1)βp−1(β − 1)k−1tp−1x(β−1)(p−2)−2

×
∣∣∣1− β + γ

β
xγ

∣∣∣p−2(
1− (β + γ)(β + γ − 1)

β(β − 1)
xγ

)]
.

If ( β(β − 1)

(β + γ)(β + γ − 1)

)1/γ

< x < 1 ,

then (
1− (β + γ)(β + γ − 1)

β(β − 1)
xγ

)
< 0

which together with (1− xγ) > 0 ensures that f(x, t) > 0. On the other hand, if

0 < x ≤
( β(β − 1)

(β + γ)(β + γ − 1)

)1/γ

,

we have (1− xγ) > k1 ≡ const. and, hence, we may find t0 small enough such that[
(1− xγ)− (p− 1)βp−1(β − 1)k−1tp−1x(β−1)(p−2)−2

×
∣∣∣1− β + γ

β
xγ

∣∣∣p−2(
1− (β + γ)(β + γ − 1)

β(β − 1)
xγ

)]
> 0 .


188 P. GIRG, L. KOTRLA, A. ŠVANDOVÁ EJDE-2022/25/CONF/26

Recall that the constants β and γ are chosen such that (β− 1)(p− 2)− 2 > 0. This
concludes the counterexample to SMP.

The counterexample above provides the existence of a nonnegative weak solution
with compact support [0, 1] for sufficiently small times. Let us observe that the solu-
tion in our counterexample is qualitatively different from the Barenblatt’s solution,
since its compact support stays constant for sufficiently small times. Moreover,
positive bump of our solution is induced by suitable nontrivial nonnegative right-
hand side of the equation rather than by the initial condition, which is identically
zero in our counterexample.

In [12], we studied counterexamples to SCP with weak solutions u and v that are
positive in (−1, 1) such that the only point satisfying u′(x) = v′(x) = 0 being x = 0.
Thus x = 0 was the only point of degeneracy of the diffusive part driven by the
p-Laplacian with p > 2. Such degeneracy causes that a perturbation from one side
of the point of degeneracy x = 0 does not spread through this point immediately in
general, but some waiting time is needed. This result from [12] is adopted to (1.1)
in the following proposition.

Proposition 3.5. Let 2 < p < ∞, max{2, p/(p− 2)} ≤ α < β, γ = β − α ( > 0),
b′(s) > k ≡ const. for all s ≥ 0, and set

u(x, t) ≡ u(x) = 1− |x|α together with

fs(x)
def
= − ∂

∂x

(∣∣∂u
∂x

∣∣p−2 ∂u

∂x

)
= (p− 1)αp−1(α− 1)|x|(α−1)(p−2)+α−2 ,

both for x ∈ (−1, 1). Furthermore, let

hi(x, t) = f(x) + |x|βψi(x, t) for (x, t) ∈ (−1, 1)× (0, T ) ; i = 1, 2,

where ψi ∈ L∞ ((−1, 1)× (0, T )) satisfy

(i) 0 ≤ ψ1 ≤ ψ2 ≤ 1
2 |x|

β(1− |x|γ + Ct) with some constant C > 0; and
(ii) ψ1 ̸≡ ψ2 on (−1, 1)× (0, τ) for any τ ∈ (0, T ).

Let wi(x, t) (i = 1, 2) be the weak solutions to

∂b(wi)

∂t
− ∂

∂x

(∣∣∂wi

∂x

∣∣p−2 ∂wi

∂x

)
= hi(x, t) for (x, t) ∈ (−1, 1)× (0, T ) ;

wi(±1, t) = 0 for t ∈ (0, T ) ;

wi(x, 0) = u0(x) for x ∈ (−1, 1)

(3.2)

which satisfy wi ∈ C ([−1, 1]× [0, T ]). If the constant C > 0 is sufficiently small,
there exists t0 > 0 such that

w1(x, t) ≤ w2(x, t) for (x, t) ∈ (−1, 1)× (0, t0) ,

w1 ̸≡ w2 on (−1, 1)× (0, t0) ,

but w1(0, t) = w2(0, t) for all t ∈ (0, t0) .

(3.3)

Finally, if ψ1 ≡ ψ2 ≡ 0 on (−1, 0) × (0, t0), then we have also w1 ≡ w2 ≡ u on
(−1, 0]× (0, t0) and w1 ̸≡ w2 on (0, 1)× (0, t0).

The proof follows the same steps as [12, Proof of Theorem 2.4, pp. 369–370].
Define

v̂(x, t) = u(x) + t |x|β (1− |x|γ) for (x, t) ∈ (−1, 1)× (0, T ),

ĝ(x, t) = b′(v̂)
∂v̂

∂t
− ∂

∂x

(∣∣∂v̂
∂x

∣∣p−2 ∂v̂

∂x

)
.


EJDE-2022/2025/CONF/26 FLOW IN POROUS MEDIA AND CFD SIMULATIONS 189

Since b′(s) ≥ k for all s ≥ 0,

ĝ(x, t)− fs(x) = b′(v̂)
∂v̂

∂t
− ∂

∂x

(∣∣∂v̂
∂x

∣∣p−2 ∂v̂

∂x

)
+

∂

∂x

(∣∣∂u
∂x

∣∣p−2 ∂u

∂x

)
≥ k

∂v̂

∂t
− ∂

∂x

(∣∣∂v̂
∂x

∣∣p−2 ∂v̂

∂x

)
+

∂

∂x

(∣∣∂u
∂x

∣∣p−2 ∂u

∂x

)
and we obtain by similar calculation as in [12, Example 2.3, pp. 368–369] that there
exists t0 ∈ (0, T ) such that

ĝ(x, t)− fs(x) ≥
1

2
|x|β

(
1− |x|γ +

C

k
t
)

for (x, t) ∈ (−1, 1)× (0, t0) ,

provided C > 0 is chosen small enough. We refer reader to [9, Proposition 1, p. 28]
for the Weak Comparison Principle (WCP) applicable to problems with b(s) ̸≡ s.
This ends the outline of the proof.

Results presented in this section demonstrate that the validity of SMP for quasi-
linear parabolic problems is a complex issue. The relevance of SMP gains impor-
tance in connection with the fact that explicit forms of solutions to problem (1.1)
can be found only in very exceptional cases and we are limited to the use of nu-
merical methods. While SMP might not hold in general, classical methods such
as implicit Euler or Crank-Nicholson methods remain reasonable effective when it
is known to be valid. However, for solutions with compact support, specialized
numerical methods are needed to mitigate unwanted numerical dispersion at the
boundary of the compact supports and avoid oscillations of the numerical solution,
see, e.g., [2, 13, 16, 23, 27, 36, 43]. Our further research on the validity of SMP
aims to provide criteria for selecting accurate and efficient numerical methods for
a given problem based on SMP validity or presence of compact support solutions.

4. Parameter estimation and CFD simulations of experiments

Understanding and managing groundwater flow in fractured hard-rock aquifers
are crucial for a variety of purposes, including sustainable water access in rural areas
and dewatering construction projects like tunnels, see, e.g., [24, 33, 50]. However,
accurately estimating parameters in constitutive laws for these aquifers presents
significant challenges compared to traditional porous media. These parameters are
usually obtained by fitting formulas such as, e.g., (2.2) or (2.3) on data obtained
by series of laboratory experiments performed on samples of given porous medium,
see, e.g., [7, 19, 28]. This process is suitable for porous media encompassing uncon-
solidated materials like soils, sands, and gravels, or permeable rocks like sandstones.
In the context of fractured hard rock, however, this approach encounters significant
limitations. Hard rocks have extremely low permeability and the groundwater flow
occurs almost exclusively within the network of fractures. It may occur that the
network of fractures can have a low density of significant fractures, making it diffi-
cult to obtain representative samples for traditional laboratory experiments. Such
experiments would require large, sometimes even multi-cubic meter rock samples to
ensure enough fractures are captured, necessitating careful extraction techniques to
minimize damage to the natural fracture network. This poses significant logistical
and cost limitations. For fracture systems, a common approach is to conduct phys-
ical experiments or numerical simulations on a single fracture or a small number


190 P. GIRG, L. KOTRLA, A. ŠVANDOVÁ EJDE-2022/25/CONF/26

of intersecting fractures, allowing the influence of fracture intersections to mani-
fest within the constitutive relationships. This strategy has been adopted in nu-
merous studies across experimental, theoretical, and numerical domains, see, e.g.,
[14, 15, 18, 32, 45, 47, 57, 60, 63].

While laboratory experiments remain valuable, approach based on computer
fluid dynamics (CFD) simulations of physical experiments offers a potentially cheaper
and more accessible alternative. By simulating the flow in several intersecting frac-
tures, we aim to estimate the parameters. To obtain as realistic outputs as possible,
we created 2D and 3D geometrical models, see Fig. 5 and 6, of several intersect-
ing fractures found on easy accessible granite rock exposures in abandoned and
partially flooded quarry Špic by Neč́ın, Czech Republic, see Fig. 3.

Figure 3. Granite rock exposures in abandoned quarry Špic by
Neč́ın, Czech Republic.

Using these geometrical models, see [58], we simulated physical experiments to
obtain datasets suitable for parameter estimation by numerical solving the station-
ary incompressible Navier-Stokes equation together with the continuity equation

−ν∆v⃗ + v⃗ · ∇v⃗ = −∇P (4.1)

∇ · v⃗ = 0 (4.2)

in the fracture network represented by the domain ΩF2D ⊂ R2 (Fig. 5) or ΩF3D ⊂
R3 (Fig. 6). Here v⃗ is the velocity field (ΩF2D → R2 or ΩF3D → R3), P is the pres-
sure field, and ν is the kinematic viscosity. The system (4.1)-(4.2) is completed by
imposing boundary conditions, see (4.3) and (4.11) for 2D and 3D case, respectively.

Let us note that our aim was to test feasibility of this approach. We do not claim
that obtained results from this preliminary stage are ready to be used in practice,
see discussion in Remark 4.1.


EJDE-2022/2025/CONF/26 FLOW IN POROUS MEDIA AND CFD SIMULATIONS 191

Figure 4. Fracture network in the granite rock exposure, quarry
Špic by Neč́ın.

Inlet →

Inlet →

→ Outlet

→ Outlet

Figure 5. 2D model of fracture network based on the most sig-
nificant fractures from the photograph on Fig. 4.

4.1. 2D model of fracture network. Our 2D model of fracture network is based
on the most significant fractures from the photograph on Fig. 4 and is represented
by domain ΩF2D ⊂ (0, 0.55)× (0, 0.5) (in meters) with Lipschitz boundary ∂ΩF2D.
This boundary ∂ΩF2D = Γinlet ∪ Γoutlet ∪ Γwall is a union of pairwise disjoint sets,
where Γinlet, Γoutlet, and Γwall represent inlet into the fracture network, outlet
from the fracture network, and fixed walls of the fracture network, respectively.
Our 2D model of fracture network has two inlets and two outlets to take into


192 P. GIRG, L. KOTRLA, A. ŠVANDOVÁ EJDE-2022/25/CONF/26

Figure 6. 3D model of fracture network. (a) Fracture network
based on the most significant fractures from the photograph on
Fig. 4. (b) 3D model with artificially added vertical fracture.

account possible mixing effects inside the network, see Fig. 5. More formally, we
have Γinlet = {(x, z) ∈ ∂ΩF2D : x = xinlet, z ∈ Iin1

∪ Iin2
} and Γoutlet = {(x, z) ∈

∂ΩF2D : x = xoutlet, z ∈ Iout1 ∪ Iout2}, where xinlet = 0, xoutlet = 0.55, Iinj
, Ioutj ⊂

(0, 0.5), j = 1, 2, are open intervals and Iin1 ∩ Iin2 = ∅, Iout1 ∩ Iout2 = ∅. Thus,
outer normal vector fields n⃗ on Γinlet and Γoutlet are constant fields, n⃗ = (−1, 0)
and n⃗ = (1, 0), respectively.

The flow in the domain ΩF2D was simulated in OpenFOAM using simpleFOAM
solver by numerically solving (4.1)-(4.2) with boundary conditions

v⃗ = (0, 0) on Γwall ,

v⃗ = (vinlet, 0) on Γinlet ,

P = 0, v satisfies condition described below on Γoutlet .

(4.3)

We impose fluxCorrectedVelocity outflow condition provided by OpenFOAM on
Γoutlet. In essence, this nonlocal condition acts computationally in the following
way. An initial estimate of the velocity at the outlet is obtained so that the change in
velocity across the outlet boundary is zero. This initial guess is then corrected based
on the calculated flux leaving the domain. This correction ensures that the flux
leaving the domain matches the expected flow based on the pressure and the internal
flow field. For detailed information, see the OpenFOAM documentation [42].

Simulation is performed for several values of vinlet, see Figure 7. In principle, this
is a numerical imitation of Darcy’s physical experiment with an adaptation for a
fractured rock. In his physical experiment, Darcy determined the total mechanical
energy loss ET using the hydraulic head h, which neglects kinetic energy but can be
easily measured in reality using piezometers. In a numerical experiment, we do not
have the possibility to measure the hydraulic head using piezometers, so we have to
choose a different approach. Just as in a physical experiment, we want to determine
the mechanical energy loss of the fluid flowing through the fracture network. In
our numerical experiment, we are modeling the actual flow in the fractures and we
work with actual velocity v⃗ not with the averaged velocity v⃗av. In our case, the
term representing kinetic energy 1/2 ρ v2 is not negligible, where v stands for the
magnitude of v⃗.


EJDE-2022/2025/CONF/26 FLOW IN POROUS MEDIA AND CFD SIMULATIONS 193

Since the velocity and pressure are dependent on location, we need to consider
average total mechanical energy per unit volume. For this purpose, we introduce

ET inlet =

∫
Γ̂inlet

ET dz∫
Γ̂inlet

dz
, (4.4)

ET outlet =

∫
Γ̂outlet

ET dz∫
Γ̂outlet

dz
, (4.5)

where

Γ̂inlet = {z ∈ R : (xinlet, z) ∈ Γinlet} = Iin1
∪ Iin2

,

Γ̂outlet = {z ∈ R : (xoutlet, z) ∈ Γoutlet} = Iout1 ∪ Iout2 .

Similar averaging approach was used in [37, p. 297] to define macroscopic pressure
and velocity.

Now, we obtain the constitutive relationship from numerical simulations as fol-
lows. We test the model for a range of values vinlet in (4.3) and record the values
of expression

△ET

△L
=
ET inlet − ET outlet

△L
.

Then we approximate dependence

△ET

△L
= flaw(uinlet)

by fitting parameters in the expression for flaw to simulated flow data. We used
the following expressions

flaw(uinlet) = αuinlet , with parameter α > 0, for Darcy’s law , (4.6)

flaw(uinlet) = β uγinlet , with parameters β, γ > 0, for power law. (4.7)

The constitutive law for specific discharge q = Ainletuinlet is then

q = Ainletf
−1
law

(△ET

△L

)
,

where Ainlet =
∫
Γ̂inlet

dz is the cross-section area of the inlet.

Now we make a connection of results of these simulations to groundwater flow in
fractured hard-rock aquifers. The water table is observed in vertical boreholes. Due
to the much larger diameter of the borehole compared to the fractures, the ground-
water flow experiences a sudden expansion upon entering. This rapid increase in
cross-sectional area significantly reduces flow velocity and so the contribution of
the term corresponding to kinetic energy to total mechanical energy per volume
in the borehole can be neglected. Thus the water would rise approximately to
level hinlet = ET inlet/(ϱg) and houtlet = ET outlet/(ϱg) above the vertical datum if
the fictive boreholes are located in the inlet and outlet of the fracture network,
respectively. This leads us to the following form of the constitutive law

q = Ainletf
−1
law

(
ϱ g

hinlet − houtlet
△L

)
. (4.8)

By fitting each law (4.6) and (4.7) to data collected from numerical simulations,
we obtained α = 17232.004 and β = 37969.871, γ = 1.852982, with significantly
lower root mean square error 17.1324 for the power-type law (4.7) compared to


194 P. GIRG, L. KOTRLA, A. ŠVANDOVÁ EJDE-2022/25/CONF/26

17232. x

37969.9 x1.85298

0.1 0.2 0.3 0.4 0.5

2000

4000

6000

8000

10000

12000

14000

∆ET /∆L

vinlet

vinlet ∆ET /∆L

0.0001 0.11832
0.0003 0.35250
0.0005 0.59079
0.0007 0.83353
0.0009 1.07690
0.001 1.19945
0.003 3.79922
0.005 6.82157
0.007 9.95944
0.009 13.71417
0.01 15.44159
0.03 69.37591
0.05 154.31617
0.07 270.63357
0.09 408.52049
0.1 493.65869
0.2 1971.17391
0.3 4057.61605
0.4 6960.98874
0.5 10507.2326

Figure 7. Table: collected data from numerical simulations per-
formed on 2D model. Graphs: Darcy’s and power type law fitted
to the data.

776.396 for Darcy’s law (4.6), see Fig. 7. This finding strongly suggests that the
power-type law provides a much more accurate description of fluid flow behavior
in the specific type of fracture networks simulated in this study. Using these fitted
values of β and γ in (4.8) together with Ainlet = 9.9987 · 10−3 (computed for our
geometrical model), and g = 9.8066 (conventional standard value for gravitational
acceleration)

q = 0.004815
(hinlet − houtlet

△L

)0.5397

. (4.9)

Assuming that the fractured hard-rock aquifer is homogeneous, isotropic (in the
sense of porosity due to fractures) and formed by a system of fractures of the same
type as in the studied section, we can use the two-dimensional differential form (2.4)
of constitutive law (4.9) for substitution in the balance equation (2.1) to obtain the


EJDE-2022/2025/CONF/26 FLOW IN POROUS MEDIA AND CFD SIMULATIONS 195

equation for the groundwater level

0.0347
∂ĥ

∂t
− 0.004815 div

(
ĥ|∇ĥ|−0.4603∇ĥ

)
= ĝ(x, y, t) , (4.10)

where we have used calculated value ϕeff = 0.0347 (the ratio of volume of the
fracture network and volume of the sample of the rock) for our 2D geometric model

of fracture network. Let us note that the term |∇ĥ|−0.4603∇ĥ is understood in the
sense that it returns zero vector for zero vector as input, cf (2.4).

Remark 4.1. Let us note that the hard-rock aquifers are often nor isotropic due
to prevailing orientation of the fractures nor homogeneous due to varying density
and aperture of the fractures, see, e.g., [29, 33, 40]. In our future research, we plan
to address these topics, especially the issue of anisotropy of nonlinear constitutive
laws.

4.2. 3D model of fracture network. Our 3D fracture network model also incor-
porates the most significant fractures from the photograph on Fig. 4. To achieve
a spatially representative network, a single artificial fracture was added and con-
nected to existing fractures within the rock mass. Although the precise location
was not based on a specific observation, it aligns with typical fracture distributions
observed in the granitic exposures of the quarry. The fracture model is then rep-
resented by domain ΩF3D ⊂ (0, 0.55) × (0, 0.5) × (0, 0.5) with Lipschitz boundary
∂ΩF3D. This boundary ∂ΩF3D = Γinlet∪Γoutlet∪Γwall is a union of pairwise disjoint
sets, where Γinlet, Γoutlet, and Γwall represent inlet into the fracture network, outlet
from the fracture network, and fixed walls of the fracture network, respectively.
Our 3D model of fracture network, see Fig. 6, is such that Γinlet ⊂ {(x, y, z) ∈
∂ΩF3D : x = xinlet} and Γoutlet ⊂ {(x, y, z) ∈ ∂ΩF3D : x = xoutlet} are sufficiently
regular, where xinlet = 0 and xoutlet = 0.55, so that outer normal vector fields n⃗ are
well defined on Γinlet and Γoutlet, respectively. Moreover, they are constant fields,
n⃗ = (−1, 0, 0) and n⃗ = (1, 0, 0) on Γinlet and Γoutlet, respectively.

The flow in the domain ΩF3D was simulated in OpenFOAM using simpleFOAM
solver by numerically solving (4.1)-(4.2) with boundary conditions

v⃗ = (0, 0, 0) on Γwall ,

v⃗ = (vinlet, 0, 0) on Γinlet ,

P = 0, fluxCorrectedVelocity condition for v⃗ on Γoutlet .

(4.11)

Simulation was performed for several values of vinlet, see Figure 8. The procedure
was analogous as in the 2D case, using averaged values of total energy over the
surfaces of inlet and outlet. By fitting each law (4.6) and (4.7), we obtained α =
7932.011 and β = 9899.873, γ = 1.888399. Again, we found that root mean square
error 12.44 for the power-type law (4.7) is significantly lower compared to 748.39
for Darcy’s law (4.6), see Fig. 8. This strongly suggests that the power-type law
is much more accurate for the specific type of fracture networks simulated in this
study.

Using these fitted values of β and γ in (4.8) together with Ainlet = 0.006 (com-
puted for our 3D geometrical model), and g as above, we find

q = 0.00597
(hinlet − houtlet

△L

)0.529549

. (4.12)


196 P. GIRG, L. KOTRLA, A. ŠVANDOVÁ EJDE-2022/25/CONF/26

7932.01 x

9899.87 x1.8884

0.2 0.4 0.6 0.8 1.0

2000

4000

6000

8000

10000

12000

∆ET /∆L

vinlet

vinlet ∆ET /∆L

0.0001 0.03358
0.0003 0.10272
0.0005 0.17364
0.0007 0.2462
0.0009 0.3228
0.001 0.36216
0.003 1.20726
0.005 2.21733
0.007 3.3564
0.009 4.60464
0.01 5.26071
0.03 23.40721
0.05 50.52472
0.07 85.79307
0.09 129.60561
0.1 154.31271
0.2 496.05733
0.3 1028.10647
0.4 1742.67221
0.5 2676.70004
0.6 3761.46226
0.7 5041.91289
0.8 6480.71039
0.9 8138.64571
1.0 9895.35284

Figure 8. Table: collected data from numerical simulations per-
formed on 3D model. Graphs: Darcy’s and power type law fitted
to the data.

Under the same assumption on the hard-rock aquifer (homogeneous, isotropic, and
formed by a system of fractures of the same type as in the studied section), we
obtain the equation for the groundwater level

0.04728
∂ĥ

∂t
− 0.00597 div

(
ĥ|∇ĥ|−0.4704∇ĥ

)
= ĝ(x, y, t) , (4.13)

where we have used calculated value ϕeff = 0.04728 (the ratio of the volume of the
fracture network and the volume of the sample of the rock) for our 3D geometric
model of fracture network.


EJDE-2022/2025/CONF/26 FLOW IN POROUS MEDIA AND CFD SIMULATIONS 197

5. P. Drábek and the p-Laplacian

Prior to embarking on the main discourse of this section, all the three authors
of this paper would like to extend their sincere congratulations to their esteemed
teacher, Pavel Drábek, on the occasion of his 70th birthday. They would like to
express their profound gratitude for his guidance and introduction to nonlinear
analysis, with a particular emphasis on the p-Laplacian, since the early stages of
their studies and scientific careers. One of the key themes of P. Drábek is the quest
for an analogue of the Fredholm alternative for nonlinear operators, particularly
for the p-Laplacian. This topic was a major focus of research at Prague School of
Nonlinear Analysis in the 1970s, where it was probably brought by J. Nečas, Ph.D.
advisor of S. Fuč́ık (who was later mentor and advisor of P. Drábek). Let us note
that Prague School of Nonlinear Analysis (or Prague School for short) was an infor-
mal research group of mathematicians primarily from Charles University Prague,
the Czechoslovak Academy of Sciences, and other academic institutions based in
Prague. The groundbreaking results of the Prague School were published in mono-
graph Spectral analysis of nonlinear operators by S. Fuč́ık, J. , J. Souček, V. Souček,
published in 1973, see [21]. After the premature death of S. Fuč́ık in 1979, this topic
gradually lost its importance in the Prague School. Among other things, because
it was a very difficult topic and after the publication of the above mentioned mono-
graph no more significant breakthroughs were achieved. Fortunately, this topic did
not completely disappear thanks to P. Drábek. His passion for this topic was ig-
nited under the mentorship of S. Fuč́ık, and he has been diligently researching it
ever since his diploma thesis. After accomplishing his C.Sc. degree (equivalent of
Ph.D.) at Czechoslovak Academy of Sciences, P. Drábek relocated to Plzeň, where
he got position in the department of mathematics at the College of Mechanical and
Electrical Engineering (precursor of today’s University of West Bohemia). Here he
continued to pursue his research in Fredholm alternative for nonlinear operators
and the p-Laplacian while the Prague School’s primary research focus shifted to
Navier-Stokes equations at the prompting of J. Nečas. P. Drábek always returned
to the topic with a certain time lag and still does. During almost half a century, an
interesting and extensive series of articles has been written, from which each rep-
resenting a major advance. In addition, he has long been attracted to this subject
attention, and so some of the seminal articles have been written without his direct
input (co-authorship), but it would hardly be without his persistence in present-
ing papers at conferences, seminars and in discussions with a number of eminent
mathematicians.

Acknowledgements. P. Girg, L. Kotrla and A. Švandová were supported by
the Grant Agency of the Czech Republic, Grant No. 22-18261S. Authors acknowl-
edge the helpful assistance of Gemini (formerly Bard), a large language model from
Google AI, in enhancing the language quality of this paper.

References

[1] Aravin, V. I.,; Numerov, S. N.; Teoriya dvizheniya zhidkostei i gazov v nedeformirue-

moi poristoi srede’. Gosudarstv. Izdat. Tehn.-Teor. Lit., Moscow, 1953. English Transl. by
A. Moscona: Theory of Fluid Flow in Undeformable Porous Media, Israel Program for Sci-
entific Translations, Jerusalem, 1965.


198 P. GIRG, L. KOTRLA, A. ŠVANDOVÁ EJDE-2022/25/CONF/26

[2] Arbogast, T.; Huang, C.-S.; Zhao, X.; Finite volume WENO schemes for nonlinear para-

bolic problems with degenerate diffusion on non-uniform meshes. Journal of Computational

Physics 399 (2019), 108921.
[3] Ayraud, V.; Aquilina, L.; Labasque, T.; Pauwels, H.; Molenat, J.; Pierson-Wickmann, A.-C.;

Durand, V.; Bour, O.;Tarits, C.; Le Corre, P.; Fourre, E., Merot, P.; Davy, P.; Compartmen-

talization of physical and chemical properties in hard-rock aquifers deduced from chemical
and groundwater age analyses. Applied Geochemistry 23, 9 (2008), 2686–2707.

[4] Barenblatt, G. I.; On self-similar motions of compressible fluid in a porous medium. Akad.

Nauk SSSR. Prikl. Mat. Meh. 16, 6 (1952), 679–698. In Russian.
[5] Barenblatt, G. I.; On some unsteady motions of a liquid and gas in a porous medium. Akad.

Nauk SSSR. Prikl. Mat. Meh. 16 (1952), 67–78. In Russian.

[6] Basak, P.; An analytical solution for the transient ditch drainage problem. Journal of Hy-
drology 41, 3 (1979), 377–382.

[7] Bear, J.; Dynamics of Fluids in Porous Media. Enviromental science series. American Elsevier
Publishing Company, Inc., New York, 1972.

[8] Bear, J.; Dynamics of Fluids in Porous Media. Dover Civil and Mechanical Engineering

Series. Dover Publications, Inc., New York, 2014.
[9] Benedikt, J.; Girg, P.; Kotrla, L.; Nonlinear models of the fluid flow in porous media and their

methods of study. In Functional differential equations and applications, FDEA-2019. Pro-

ceedings of the 7th international conference, Ariel, Israel, September 22–27, 2019. Singapore:
Springer, 2021, pp. 15–42.

[10] Benedikt, J.; Girg, P.; Kotrla, L.; Takáč, P.; The strong maximum principle in parabolic

problems with the p-Laplacian in a domain. Appl. Math. Lett. 63 (2017), 95–101.
[11] Benedikt, J.; Girg, P.,;Kotrla, L.; Takáč, P.; Origin of the p-Laplacian and A. Missbach.

Electron. J. Differential Equations (2018), Paper No. 16, 17 pp.

[12] Benedikt, J.; Girg, P.; Kotrla, L.; Takáč, P.; The strong comparison principle in parabolic
problems with the p-Laplacian in a domain. Appl. Math. Lett. 98 (2019), 365–373.

[13] Berger, Alan E.; Brezis, H.; Rogers, J. C. W.; A numerical method for solving the problem
ut − δf(u) = 0. ESAIM: Mathematical Modelling and Numerical Analysis - Modélisation

Mathématique et Analyse Numérique 13, 4 (1979), 297–312.

[14] Berkowitz, B.; Characterizing flow and transport in fractured geological media: A review.
Advances in water resources 25, 8-12 (2002), 861–884.

[15] Brush, D. J.; Thomson, N. R.; Fluid flow in synthetic rough-walled fractures: Navier-stokes,

stokes, and local cubic law simulations. Water Resources Research 39, 4 (2003).
[16] Cavalli, F., Naldi, G.; Puppo, G.; Semplice, M.; High-order relaxation schemes for nonlinear

degenerate diffusion problems. SIAM Journal on Numerical Analysis 45, 5 (2007), 2098 –

2119.
[17] Cheng, H.; Wang, F.; Li, S.; Guan, X.; Yang, G.; Cheng, Z.; Yu, C.; Yuan, Y.; Feng, G.;

Effect of movability of water on the low-velocity pre-darcy flow in clay soil. Journal of Rock

Mechanics and Geotechnical Engineering (2024).
[18] Cherubini, C.; Giasi, C.; Pastore, N.; Evidence of non-darcy flow and non-fickian transport

in fractured media at laboratory scale. Hydrology and Earth System Sciences 17, 7 (2013),
2599–2611.

[19] Darcy, H.; Les fontaines publiques de la ville de Dijon. Victor Dalmont, Paris, 1856.

[20] Forchheimer, P.; Wasserbewegung durch Boden. Zeit. Ver. Deutsch. Ing. 45 (1901), 1736–
1741 and 1781–1788.

[21] Fuč́ık, S.; Nečas, J.; Souček, J.; Souček, V.; Spectral analysis of nonlinear operators. Lecture
Notes in Mathematics, Vol. 346. Springer-Verlag, Berlin-New York, 1973.

[22] Geoscience Australia; http://www.ga.gov.au/scientific-topics/water/groundwater/

groundwater-in-australia/fractured-rocks, 2014. [Online; accessed 21-February-2020].

[23] Gu, Y.; Shen, J.; Bound preserving and energy dissipative schemes for porous medium equa-
tion. Journal of Computational Physics 410 (2020), 109378.

[24] Gustafson, G.; Krásný, J.; Crystalline rock aquifers: Their occurrence, use and importance.
Applied Hydrogeology 2, 2 (1994), 64–75.

[25] Harr, M.; Groundwater and Seepage. Dover Civil and Mechanical Engineering. Dover Publi-

cations, 2012.
[26] Izbash, S. V.; O filtracii v krupnozernistom materiale. Izv. Nauchno-Issled. Inst. Gidro-Tekh.

(N.I.LG.), Leningrad 1, 1931 (In Russian).

http://www.ga.gov.au/scientific-topics/water/groundwater/groundwater-in-australia/fractured-rocks
http://www.ga.gov.au/scientific-topics/water/groundwater/groundwater-in-australia/fractured-rocks


EJDE-2022/2025/CONF/26 FLOW IN POROUS MEDIA AND CFD SIMULATIONS 199

[27] Jiang, Y.; High order finite difference multi-resolution WENO method for nonlinear degen-

erate parabolic equations. Journal of Scientific Computing 86, 1 (2021).

[28] King, F.; Principles and conditions of the movements of ground water. Nineteenth Ann. Kept.
U. S. Geol. Survey pt. 2, 9-12 (1898), 209–215.

[29] Kiràly, L.; Groundwater flow in heterogeneous, anisotropic fractured media: A simple two-

dimensional electric analog. Journal of Hydrology 12, 3 (1971), 255–261.
[30] Kröber, C.; Versuche über die bewegung des wassers durch sandschichten. Zeitschr. des Vere-

ines deutscher Ing. 28, 31 and 32 (1884), 593–595 and 617 –619.

[31] Lachassagne, P.; Dewandel, B.; Wyns, R.; Review: Hydrogeology of weathered
crystalline/hard-rock aquifers—guidelines for the operational survey and management of their

groundwater resources. Hydrogeology Journal 29, 8 (Dec 2021), 2561–2594.

[32] Lapcevic, P. A.; Novakowski, K. S.; Sudicky, E. A.; The interpretation of a tracer experiment
conducted in a single fracture under conditions of natural groundwater flow. Water Resources

Research 35, 8 (1999), 2301–2312.
[33] Larsson, I.; Unesco; 8.6, I. H. P. P.; Ground Water in Hard Rocks: Project 8.6 of the

International Hydrological Programme. Studies and reports in hydrology. Unesco, 1984.

[34] Leibenson, L. S.; Turbulent movement of gas in a porous medium. Bull. Acad. Sci. USSR.
Sér. Géograph. Géophys. [Izvestia Akad. Nauk SSSR] 9 (1945), 3–6. In Russian. Reprinted

in Ref. [35], 499–502.

[35] Leibenson, L. S.; Sobranie trudov, Chast’ II: Podzemnaya gidrodinamika [Collected Works,
Vol. II: Underground Hydrodynamics]. Izdat’elstvo Akademii Nauk S.S.S.R., Moscow,

U.S.S.R., 1953. In Russian.

[36] Liu, Y.; Shu, C.-W.; Zhang, M.; High order finite difference WENO schemes for nonlinear
degenerate parabolic equations. SIAM Journal on Scientific Computing 33, 2 (2011), 939 –

965.

[37] Lucas, Y.; Panfilov, M.; Buès, M.; High velocity flow through fractured and porous media:
the role of flow non-periodicity. European Journal of Mechanics - B/Fluids 26, 2 (2007),

295–303.
[38] Macdonald, A.; Davies, J.; Calow, R.; African hydrogeology and rural water supply. Applied

Groundwater Studies in Africa (2008), 127–148.

[39] Marino, M.; Rise and decline of the water table induced by vertical recharge. Journal of
Hydrology 23, 3-4 (1974), 289–298.

[40] Maréchal, J. C.; Dewandel, B.; Subrahmanyam, K.; Use of hydraulic tests at different scales

to characterize fracture network properties in the weathered-fractured layer of a hard rock
aquifer. Water Resources Research 40, 11 (2004).

[41] Missbach, A. A.; Filtrovatelnost čeřených a saturovaných št’áv. IV. Přezkoušeńı vzorce van

Gilse, ... Listy cukrov. 54, 39 (1936), 361 – 368. In Czech.
[42] Open Foam; https://www.openfoam.com. [Online; accessed 02/28/2024 23:36].

[43] Parlange, J.-Y.; Hogarth, W.; Govindaraju, R.; Parlange, M.; Lockington, D.; On an exact

analytical solution of the Boussinesq equation. Transport in Porous Media 39, 3 (2000), 339
– 345.

[44] Perrin, J.; Ahmed, S.; Hunkeler, D.; The effects of geological heterogeneities and piezomet-
ric fluctuations on groundwater flow and chemistry in a hard-rock aquifer, southern India.

Hydrogeology Journal 19, 6 (2011), 1189–1201.

[45] Quinn, P.; Cherry, J.; Parker, B.; Relationship between the critical Reynolds number and
aperture for flow through single fractures: Evidence from published laboratory studies. Jour-

nal of Hydrology 581 (2020), 124384.
[46] Rai, S.; Singh, R.; Two-dimensional modelling of water table fluctuation in response to

localised transient recharge. Journal of Hydrology 167, 1 (1995), 167–174.

[47] Sarkar, S.; Toksoz, M. N.; Burns, D. R.; Fluid flow modeling in fractures. Tech. rep., Mas-

sachusetts Institute of Technology. Earth Resources Laboratory, 2004.
[48] Scheidegger, A. E.; The Physics of Flow through Porous Media. The Macmillan company,

New York, 1960.
[49] Şen, Z.; Applied Hydrogeology for Scientists and Engineers. CRC Press, Boca Raton, 1995.
[50] Shapiro, A. M.; Fractured-rock aquifers understanding an increasingly important source of

water. https://pubs.usgs.gov/publication/fs11202, 2002. [Online; accessed 02/28/2024

23:43].

 https://www.openfoam.com
https://pubs.usgs.gov/publication/fs11202


200 P. GIRG, L. KOTRLA, A. ŠVANDOVÁ EJDE-2022/25/CONF/26

[51] Singh, R.; Rai, S.; On subsurface drainage of transient recharge. Journal of Hydrology 48,

3-4 (1980), 303–311.

[52] Singh, R.; Rai, S.; A solution of the nonlinear Boussinesq equation for phreatic flow using an
integral balance approach. Journal of Hydrology 109, 3-4 (1989), 319–323.

[53] Singhal, B. B. S.; Nature of Hard Rock Aquifers: Hydrogeological Uncertainties and Ambi-

guities. Springer Netherlands, Dordrecht, 2008, pp. 20–39.
[54] Smreker, O.; Entwicklung eines Gesetzes für den Widerstand bei derBewegung des Grund-

wassers. Zeitschr. des Vereines deutscher Ing. 22, 4 and 5 (1878), 117–128 and 193–204.

[55] Smreker, O.; Das grundwasser und seine verwendung zu wasserversorgungen. Zeitschr. des
Vereines deutscher Ing. 23, 4 (1879), 347–362.

[56] Soni, J.; Islam, N.; Basak, P.; An experimental evaluation of non-darcian flow in porous

media. Journal of Hydrology 38, 3–4 (1978), 231–241.
[57] Stark, K. P.; Volker, R. E.; A study of some theoretical aspects of non-linear flow through

porous materials. Tech. Rep. in Res. Bulletin No. 1, Dept. of Civil Engineering, University
College of Townsville, Townsville, Australia, April 1967.

[58] Švandová, A.; Simulation of fluid flow in porous media with an emphasis on fractured me-

dia and coarse-grained materials. Master’s thesis, University of West Bohemia, 2022. Men-
tor: P. Girg.

[59] Wright, E.; The hydrogeology of crystalline basement aquifers in Africa. Geological Society

Special Publication 66 (1992), 1–27.
[60] Yan, X.; Qian, J.; Ma, L.; Wang, M.; Hu, A.; Non-fickian solute transport in a single fracture

of marble parallel plate. Geofluids 2018 (2018).

[61] Zhukovskii, N. E.; Teoreticheskoe issledovanie o dvizhenii podpochvennykh vod. Zhurnal
Russkogo fiziko-khimicheskogo obshchestva 21, 1 (1889). [In Russian], reprinted in Ref. [62],

p. 9 – 33.

[62] Zhukovskii, N. E.; Polnoe sobranie sochinenii, vol. 7. Moscow-Leningrad, 1937. [Collected
Papers], Russian with English Summary.

[63] Zimmerman, R. W.; Al-Yaarubi, A.; Pain, C. C.; Grattoni, C. A.; Non-linear regimes of fluid
flow in rock fractures. International Journal of Rock Mechanics and Mining Sciences 41,

SUPPL. 1 (2004), 163–169.

[64] Zunker, F.; Das allgemeine Grundwasserfliessgesetz. Journal für Gasbeleuchtung und
Wasserversorgung 63, 21 (1920), 331–334, and 350.

Petr Girg

Department of Mathematics and NTIS, Faculty of Applied Scences, University of West

Bohemia, Univerzitńı 8, CZ-301 00 Plzeň, Czech Republic
Email address: pgirg@kma.zcu.cz

Lukáš Kotrla
Department of Mathematics and NTIS, Faculty of Applied Scences, University of West

Bohemia, Univerzitńı 8, CZ-301 00 Plzeň, Czech Republic

Email address: kotrla@ntis.zcu.cz

Anežka Švandová
Department of Mathematics, Faculty of Applied Scences, University of West Bohemia,
Univerzitńı 8, CZ-301 00 Plzeň, Czech Republic

Email address: svandova@kma.zcu.cz


	1. Introduction
	2. Mathematical models
	3. Maximum and comparison principles
	4. Parameter estimation and CFD simulations of experiments
	4.1. 2D model of fracture network
	4.2. 3D model of fracture network.

	5. P. Drábek and the p-Laplacian
	Acknowledgements

	References