context of porous media where $\Delta s$ is the Laplacian operator applied to $s$ and $T$ represents the effective surface tension. This term is analogous to the linearized Laplace pressure of a thin liquid film. Recently, [20] proposed a more complex expression for the surface energy. In our context, such a surface energy term is crucial since it plays a regularizing role. In the absence of this effective surface tension $T$, our model with hydrophobicity leads to infinitely small structures and a diverging growth rate as shown later (see Eq. (10)). In this paper we focus on the simplest model in [15] to reproduce fingering and water trapping. Therefore, the flux $\vec{\Phi}$ is now given by:

$$ \vec{\Phi} = k_r(s) \nabla (\Pi(s) + z + T\Delta s). \quad (6) $$

In the case of a stratified medium, we assume that the properties depend only on the depth $z$. In the hydrophilic layer, the medium follows the Genuchten-Mualem law, i.e. $h_f = 0$ in our model Eq. 4, and in the hydrophobic layer $h_f = h_f^o \neq 0$. We introduce the function $\xi(x, z)$ that realizes a continuous and derivable transition of the parameter $h_f$ between the different layers such that $h_f(x, z) = \xi(x, z)h_f^o$ and $\xi(x, z) \simeq 0$ in the hydrophilic layer while $\xi(x, z) \simeq 1$ in the hydrophobic layer. Then the wettability function reads

$$ \Pi(s, x, z) = h_c(s) - \xi(x, z)h_r(s). \quad (7) $$

The equation can be written:

$$ \frac{\partial s}{\partial t} = \nabla \cdot \left[ k_r(s) \nabla \left( \frac{\delta \mathcal{F}}{\delta s} \right) \right] - \frac{\partial k_r(s)}{\partial z} \quad (8) $$

where $\frac{\delta \mathcal{F}}{\delta s} = -\Pi(s, z) - T\Delta s$ is the functional derivative of the functional $\mathcal{F}$:

$$ \mathcal{F}(s) = \int_D \left( \left( \int_s^1 \Pi(s, z) ds \right) - T \frac{1}{2} |\nabla s|^2 \right) dx dz. \quad (9) $$

The functional $\mathcal{F}$ is also a Lyapunov functional and represents the free energy of the system [21]. Eq. (8) is a fourth order and nonlinear Partial Differential Equation (PDE). This equation is reminiscent of the lubrication equation with wettability in thin film theory [14]. We then expect front instability and intermittent flow as shown in [22,23].

### Horizontal diffusion in a hydrophobic medium.

We consider the 1D model of a horizontal imbibition with periodic boundary conditions at $x = \pm L/2$ where $L$ is the domain size. For profiles that are symmetric by the reflection $x \leftrightarrow -x$, the latter condition also implies the no-flux condition $\frac{\partial s}{\partial n} = 0$ at these boundaries. In this geometry the water content is conserved and we denote $s_m$ the mean water content.

**Homogeneous medium.** The constant distribution $s \equiv s_m$ is the trivial solution of the 1D horizontal problem but this solution may be unstable because of the repellent term $h_r$. The linear stability is given by the eigenmodes of the

<figure><img src="image_42.png" /><figcaption>Fig. 2: Horizontal saturation $s(x)$ profile in hydrophobic medium for different values of mean saturation $s_m$. Plain [dashed] line indicates stable [unstable] saturation distribution. Parameters are $s_{rg} = 0.2$; $h_f = 0.9$; $\delta = 0.1$, $T = 0.001$.</figcaption></figure>

linearization of the right hand expression in (8). Because of the periodic condition and the translation invariance in the $x$ direction, the eigenmodes are the Fourier modes. A standard derivation shows that the eigenvalue of the $k$ wavenumber is $\beta(k) = k_r(s_m)k^2(\Pi'(s_m) - Tk^2)$. The film is unstable if $\beta(k) > 0$ which requires i) $\Pi'(s_m) = \Pi'_m > 0$ and ii): $k < k_c = \sqrt{\frac{\Pi'_m}{T}}$. The condition i) is realized if $s_m$ is in the interval $[s_1, s_2]$ and the condition ii) is realized for longwaves $\lambda$ such that

$$ \lambda > \lambda_c = \frac{2\pi}{k_c} = 2\pi \sqrt{\frac{T}{\Pi'_m}}. \quad (10) $$

The maximum growth rate corresponds to the wavelength $\lambda_m = \lambda_c \sqrt{2}$ [14]. When the linear instability evolves an alternation of dry and wet regions periodically appears with the characteristic length of about $\lambda_m$. Note that $T > 0$ ensures that $\lambda_m > 0$ and hence a finite size of the structure as mentioned in the description of the Modelling. The development of instability is analogous to dewetting on a substrate or spinodal decomposition [24]. The resulting dry and wet regions coalesce through a coarsening phenomenon due to the action of the macroscopic surface tension. The coalescence finishes when a unique wetted region coexists with a dry region. This steady-state corresponds to the minimum of the free energy $\mathcal{F}(s)$. We expect, then, the imbibition dynamics to converge towards this steady-state $s_{min}(x)$ as in [22]. The profiles of $s_{min}$ displayed in Fig. 2 show that the moist area is localized in a finite domain if $s_m$ is below a critical value $s_c = 0.739$. In order to better understand the origin of the transition at $s = s_c$, we analyse the bifurcation diagram of steady-state Fig. 3 by increasing the mean saturation $s_m$.

From the trivial solution $s \equiv s_m$, a solution branch emerges subcritically at $s_1 = 0.0449$. The solution is unstable till a turning point $s_{sn1} = 0.02$. The branch from this turning point corresponds to the minimum free energy and therefore it is stable. For a water content, three equilibria coexist in the narrow range $[s_1, s_{sn1}]$: The trivial solution (stable) and the pair of stable/unstable nonuniform solutions. Fig. 2 displays both nonuniform profiles for $s_m = 0.03$. The water spreads over the entire domain
<figure><img src="image_43.png" /><figcaption>Fig. 3: Horizontal saturation profile in a heterogeneous hydrophobic medium for different values of mean saturation $s_m$. Bold [fine] lines indicate stable [unstable] steady-state solutions. The horizontal line corresponds to solution branches with uniform distribution. Red symbols '+' indicate the steady-states whose profiles are depicted in Fig. 4 for $s_m = 0.03, 0.10, 0.40, 0.70$. Parameters as in Fig. 2.</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>Fig. 4: [black dashed line] Heterogenous function $\xi(x)/10$ and [color plain lines] horizontal saturation profile in an inhomogeneous medium for different values of mean saturation $s_m$. The range where $\xi(x) \simeq 0$ corresponds to the hydrophobic region while $\xi(x) \simeq 1$ is the hydrophilic region. Parameters as in Fig. 2.</figcaption></figure>

for the unstable equilibrium while for the stable equilibrium the main water content is concentrated in a small range and the saturation maximum is about 0.53.

Increasing the mean saturation along this stable branch, the water remains localized and the saturation maximum increases till $s_w = 0.78$ is reached for $s_m$ about 0.10 (Fig. 2).

For a mean saturation greater than 0.10, this time it is the saturation maximum which remains almost constant at $s_w$ and only the width of the wetted domain increases with the water content. Thus, water distribution is divided into two zones, one wet with a constant saturation about $s_w = 0.78$, the other dry with a low saturation $s_d = 0.009$ (see the profiles for $s_m = 0.40$ and $s_m = 0.70$). These saturation values are such that $\Pi(s_w) \simeq \Pi(s_d)$ which is the equilibrium condition when $T\Delta s$ is negligible. The width of the wet area increases with mean water content $s_m$ till the dry region vanishes due to the finite size domain. Then, the solution branch ends at a second turning point corresponding to $s_c = 0.739$ (Fig. 3). The turning unstable branch connects the trivial solution at $s_2 \simeq 0.47$. We represent the equilibrium for $s_m = 0.70$ of this unstable branch in Fig. 2. There is no longer a dry area and the saturation contrast between the minimum, 0.08, and the maximum, 0.73, of saturation is weaker than for the stable solution. For $s_m > s_2$ the uniform water distribution is also stable and then in the interval $[s_2, s_c]$, the trivial uniform solution coexists with the non-uniform profile and both are stable. The final state depends on the initial conditions.

Beyond $s_c$ the final state corresponds to the uniform distribution as for the hydrophilic medium.

To conclude, for a domain that is large enough the water spreads over a finite area in a homogeneous medium. In the wet area, the water saturation $s_w$ is almost constant and its value depends on the $\Pi(s)$ function. We find numerically that $s_w \simeq 0.78$ for our parameters.

Note that there are many solution branches emerging from

the points $s_1$ and $s_2$ which correspond to patterns with the periodicity $L/n$, $n \ge 2$. Nevertheless, as discussed before all these branches remain unstable.

*Inhomogeneous medium.* We consider an inhomogeneous medium consisting of a hydrophobic medium containing a hydrophilic region, i.e. $h_f = 0$. The wettability $\Pi(s, x)$ is given by Eq. 7 with

$$\xi(x) = \frac{1}{2} - \frac{1}{2} \tanh \left( 5 \cos \left( \pi \tanh \left( \frac{2\pi}{L} 1.8x \right) \right) \right) \quad (11)$$

The function $\xi(x)/10$ is displayed Fig. 4. The region where $\xi(x) \simeq 1$ corresponds to a hydrophilic medium. The steady-states for different water contents $s_m$ show that water is confined in the hydrophilic region for small $s_m$ (Fig. 4). When the maximum saturation in the hydrophilic region is about 0.50 (see the profile for $s_m = 0.13$ in Fig. 4) a small amount of water invades the hydrophobic region. Then, the water content in the hydrophobic region increases progressively with $s_m$. The saturation maximum ceases to increase when the mean water content $s_m$ reaches 0.40. Then, the saturation maximum is about 0.92 and a further increase in water content implies a spread in the hydrophobic region (Fig. 4). As for the homogeneous medium, the spread is finite and the saturation maximum is about 0.80.

The model reproduces a finite diffusion of water in homogeneous or inhomogeneous media. In the following, we focus on a vertical stratified medium.

**Stratified medium.** — In this section, we consider a periodic vertical stratification of a hydrophobic/hydrophilic medium such that the heterogeneity depends only on $z$, i.e. in the horizontal $x$ direction, the medium is assumed to be homogeneous. We focus on a spatially periodic flow in the $z$ direction. The vertical domain length is fixed to $L_z = 1$ and we apply periodic boundary conditions at $z = 0$ and $z = 1$. The function $\xi(z)$ characterizes the wettability variation and it is given by Eq. 11 replacing the variable $x$ by $z + 0.35$.
<figure><img src="image_45.png" /><figcaption>Fig. 5: Saturation profiles of steady-states. [Upper panel] Profiles from $s_m = 0.04$ to 0.185 and [lower panel] profiles from $s_m = 0.195$ to 0.60. Bold [fine] lines indicate stable [unstable] steady-states. Dashed bold lines indicate the profile at the bifurcations: [upper panel] saddle-node $s_c \approx 0.19$ and [lower panel] Hopf bifurcation $s_{Hopf} = 0.44$. The black dotted line corresponds to $\xi(x)/10$. Parameters as in Fig. 2.</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>Fig. 6: Flux $\Phi$ depending on $s_m$ for each 1D equilibrium in a vertical stratified medium. Bold plain [dashed] lines indicate stable [unstable] steady-state. The fine plain straight line is the flux associated to $k_r(s_m)$. Black dots correspond, by increasing $s_m$, to the $s_{2d}$, $s_{sn}$ and $s_h$ bifurcations. Steady-states corresponding to the profiles shown in Fig. 5 are indicated by + symbols. Mean flux ratios for time periodic flows are indicated by red dots. Parameters as in Fig. 2.</figcaption></figure>

Then the hydrophilic layer is localized in the vertical range $z \in [0, 0.3]$. In contrast to the horizontal case, the dynamics is no longer a relaxation dynamics to the minimum of $F$ but we expect front instability and intermittent flow as shown in [23].

**1D vertical stratified medium.** Because of the periodic conditions, the mean water content noted $s_m$ is constant. We study steady-state solutions by increasing the water content $s_m$. Compared to the horizontal case, the maximum of water saturation is shifted downstream near the interface between the hydrophilic/hydrophobic layers. For low mean saturation, the hydrophilic layer can retain most of the water and the water is blocked at the interface resulting in a stationary state displayed in Fig. 5 (upper panel). On increasing the mean saturation $s_m$, the saturation maximum approaches the interface and a small amount of water spreads downwards into the hydrophobic region.

The bifurcation diagram Fig. 6 represents the flux $\Phi$ of the latter steady-states as a function of $s_m$. The solution branch displays a first saddle-node at $s_m = s_{sn} = 0.197383$. For a saturation near the saddle-node $s_m > 0.19$, the saturation profile of the steady-state changes noticeably. The saturation maximum decreases and a greater amount of water advances significantly into the hydrophobic part. The latter steady-states are stable till the value $s_h = 0.197345$ corresponding to a Hopf bifurcation. This bifurcation leads to a subcritical branch of unstable small amplitude time periodic states. We argue that the bifurcation scenario is similar to the one described in [23] on ridge depinning on a wettability defect. The unstable branch of periodic orbit undergoes a global bifurcation at

$s_g = 0.197170 < s_h < s_{sn}$ (not shown in Fig. 6) leading to an intermittent flow.

The time integration for $s_m > s_g$ shows a time-dependent dynamics. When $s_m$ is close to the global bifurcation, there are two main steps inducing an intermittent-like dynamics. In the first phase, the water from the hydrophilic layer spreads slowly in the hydrophobic region. There are two local maxima of the saturation profile, one near the interface and one at the head of the front in the hydrophobic region ($t = 667$ in Fig. 7). Between the two maxima the saturation is lower than $s_2$ implying a hydrophobic behavior of the medium. This results in a decrease in the water content in this region followed by the separation of the downstream water in the hydrophobic layer from the water remaining in the hydrophilic layer ($t = 671$). The water in the hydrophobic layer advances like a solitary wave and in its wake leaves a small amount of water, about 5% ($t = 671-674$). This wave joins the next hydrophilic layer and is rapidly absorbed by the hydrophilic layer where the water remains trapped for a considerable length of time: $\Delta t \approx 500$. Thus, the water infiltrates the hydrophobic layer in a similar manner and an intermittent-like dynamics takes place. Note that for mean saturation $s_m$ in the range [$s_g, s_{h2}$] there also exists a steady-state flow represented in Fig. 6 by the dashed line. However, this stationary flow is linearly unstable and the asymptotic dynamics remains the intermittent dynamics.

On increasing the mean saturation $s_m$, the dynamics loses its intermittent character and instead periodic oscillations of water content take place between the two layers. This branch of solutions ends at a Hopf bifurcation ($s_{h2} = 0.44$) with the steady-state solution branch. The profile resembles a wavy profile with higher water content in the hy-
<figure><img src="image_47.png" /><figcaption>Fig. 7: Profiles of water saturation at different times *t* during the intermittent flow for *s*<sub>m</sub> = 0.1972. Parameters as in Fig. 2.</figcaption></figure>

drophilic layer but both layers are moist (see the bottom panel in Fig. 5). The steady-state solution is stable for *s*<sub>m</sub> > *s*<sub>h2</sub> and remains stable till complete saturation of the medium.

**2D geometry: Fingering instability.** As for the instability analysis in a homogeneous medium, a longwave instability may arise in the transverse direction *x* if the lateral length *ℓ* is large enough. Since the problem is invariant in the *x* direction, the eigenmodes take the form *u*(*x*, *z*, *t*) = *u*<sub>*x*</sub>(*z*) exp(−*ik*<sub>*x*</sub>*x* + *β*(*k*<sub>*x*</sub>)*t*), where *k*<sub>*x*</sub> is the wave number in the *x*-direction and *β*(*k*<sub>*x*</sub>) the growth rate. In the limit *k*<sub>*x*</sub> → 0, every function *u*<sub>*x*</sub>(*z*) is an eigenmode from the 1D case except for one eigenmode called varicose and such that *u*<sub>*x*</sub>(*z*) is not zero averaged over *z* [25]. Since the mass is conserved in our problem, the eigenvalue associated with this mode is zero for *k*<sub>*x*</sub> = 0. Thus, we study the transverse stability from a stable 1D stationary state. The varicose mode may become strictly positive when *k*<sub>*x*</sub> becomes strictly positive, i.e. it is at the origin of the longwave instability [14]. Fig. 8 shows the eigenvalue *β*(*k*<sub>*x*</sub>) of the varicose mode for two steady-states from the bifurcation diagram 6 at *s*<sub>*m*</sub> = 0.16 and *s*<sub>*m*</sub> = 0.1798. For *s*<sub>*m*</sub> = 0.16, the eigenvalues become strictly negative when *k*<sub>*x*</sub> > 0 and then the steady-state is also stable for the 2D problem. In contrast for *s*<sub>*m*</sub> = 0.1798, the eigenvalue *β*(*k*<sub>*x*</sub>) is strictly positive in the interval [0, *k*<sub>*c*</sub>] (*k*<sub>*c*</sub> = 7.24) and reaches the maximum value *β*<sub>*max*</sub> = 0.0967 at *k*<sub>*max*</sub> = 4.1 and *ℓ*<sub>*max*</sub> = 1.53.

The lateral instability of the 1D-steady-state is shown in
red in the bifurcation diagram 6 corresponding to the in-
terval [$s_{\ell}$; $s_{sn}$]. This instability leads to the spontaneous
formation of fingers as shown by the time integration start-
ing from the 1D profile for $s_m = 0.1798$. Adding a sinu-
soidal perturbation of 1% amplitude with period $\ell = 2$,
we observe the protrusion growing in a finger (Fig. 9) at
t = 15. This time is consistent with the characteristic
time $t = 1/\beta = t_c \approx 12$ of the linear growth. The water
at the base of the protrusion recedes to the hydrophilic
layer while the finger stretches ($t = 21$ in Fig. 9). At the
fingerhead, we observe a oversaturation which corrobo-
rates the experimental observations in porous media [26].

<figure><img src="image_48.png" /><figcaption>Fig. 8: Dispersion relation of the 1D steady-state in stratified medium for $s_m = 0.16$ (fine line) and $s_m = 0.1798$ (bold line). Parameters as in Fig. 2.</figcaption></figure>

In the tail of the finger the width decreases due to the
hydrophilicity (t = 23). Thus a small amount of water
like a drop is almost detached from the hydrophilic layer.
However, the humidity in its wake is not negligible, cre-
ating a preferential path. This water travels through the
hydrophobic layer and reaches the next hydrophilic layer
(t = 24). Before the water has been completely absorbed
by the hydrophilic part, a new finger appears in the previ-
ous finger path (t = 24.5). This scenario repeats but after
each cycle the water saturation in the path increases until
it converges to a stationary state (t = 61.75): most of the
water is in the hydrophilic layer while a narrow rivulet
connects these layers.

This time simulation is reminiscent of experiments dis-
playing finger flows in a hydrophobic sand [11, 27] not
only qualitatively (see, for instance, Fig. 6 in [11]) but
also quantitatively. The finger width *w* ranges from 1 to
5 cm. Considering a grain size of about 0.3 *mm*, the capil-
lary height *H*<sub>c</sub> is about 20 *cm* according to [28]. Then the
relative finger width *w*/*H*<sub>c</sub> ranges from 0.05 to 0.25 which
fits with the width 0.2 of the simulated rivulet.

**Concluding remarks.** — Taking into account sur-
face energies especially wettability by a disjoining pres-
sure in a variant model of the Darcy-Richards equation, we
capture water trapping in hydrophilic layers, intermittent-
like flow and fingering phenomena in a stratified medium.
The keypoint of the modeling, by adding a hydrophobic
term, is that the suction function is no longer monotonous.
The hydrophobic term represents the effective action of
hydrophobic molecules and then depends on the satura-
tion. This term is added to the standard capillary pressure
and there is a range over which hydrophobicity is 'visible'
which corresponds to the positive slope of Π(s). According
to [15], we add a term related to the surface which takes
water free surface energy into account. In a more general
way, Eq. 8 with the free energy term (Eq. 9) belongs to
mesoscopic models of complex wetting [21].

Using bifurcation analysis tools, we show that the 1D flow
threshold results from a global bifurcation. This analysis
has not only a theoretical interest but it makes it possible
to determine the logarithmic time scale governing inter-
Model of hydrophobic stratified porous media

<figure><img src="image_49.png" /><figcaption>Fig. 9: (Colour on-line) Snapshots of the time evolution of water infiltration from hydrophobic layer into the hydrophilic layer. The initial condition is the 1D steady-state for $s_m = 0.1798$ (see Fig. 6) perturbed by the sinusoidal function $\cos(\pi x)$ with a relative amplitude of 1%. Colour code of the water content is identical to all the panels. Parameters as in Fig. 2.</figcaption></figure>

mittency. Moreover, the spontaneous formation of finger-
ing constitutes a preferential path for water, inducing a
hysteresis effect. This bifurcation scenario is similar to
the one described in depth in [23] in the thin film flow
context.

As long as the wetting is only partial, the current model is
relevant. Therefore, the instability term introduced in [20]
to simulate fingering in sand can also be interpreted as a
partial wetting property of sand.

REFERENCES

[1] DOERR S., SHAKESBY R. and WALSH R., *Earth-Science Reviews*, **51** (2000) 33.

[2] DICARLO D. A., *Water Resources Research*, **49** (2013) 4531.

[3] HILL D. E. and PARLANGE J.-Y., *Soil Sci. Soc. Am. J.*, **36** (1972) 697-702.

[4] PARLANGE J.-Y. and HILL D. E., *Soil Science*, **122** (1976).

[5] HASSANIZADEH S. M. and GRAY W. G., *Water Resources Research*, **29** (1993) 3389.

[6] ELIASSI M. and GLASS R. J., *Water Resources Research*, **38** (2002) 16.

[7] CUETO-FELGUEROSO L. and JUANES R., *Phys. Rev. E*, **79** (2009) 036301.

[8] RÄTZ A. and SCHWEIZER B., *ZAMM - Journal of Applied Mathematics and Mechanics*, (2012).

[9] DEKKER L. W. and RITSEMA C. J., *Water Resources Research*, **30** (1994) 2507.

[10] RITSEMA C. J. and DEKKER L. W., *Water Resources Research*, **30** (1994) 2519.

[11] NIEBER J., BAUTERS T., STEENHUIS T. and PARLANGE J.-Y., *Journal of Hydrology*, **231-232** (2000) 295.

[12] BELIAEV A. and HASSANIZADEH S., *Transport in Porous Media*, **43** (2001) 487.

[13] LIU H., JU Z., BACHMANN J., HORTON R. and REN T., *Soil Sci. Soc. Am. J.*, **76** (2012) 342.

[14] KALLIADASIS S. and THIELE U., (Editors) *Thin films of soft matter* (Springer, Wien) 2007 Ch. Structure formation in thin liquid films pp. 25–94.

[15] CUETO-FELGUEROSO L. and JUANES R., *Water Resources Research*, **45** (2009) W10409.

[16] VAN GENUCHTEN M. T., *Soil Science Society American Journal*, **44** (1980) 892.

[17] VEREECKEN H., WEYNANTS M., JAVAUX M., PACHEP-SKY Y., SCHAAP M. G. and VAN GENUCHTEN M., *Vadose Zone Journal*, **9** (2010) 795.

[18] THIELE U., VELARDE M. G., NEUFFER K. and POMEAU Y., *Phys. Rev. E*, **64** (2001) 031602.

[19] BELTRAME P., *Europhys. Lett.*, **121** (2018) 64002.

[20] BELJADID A., CUETO-FELGUEROSO L. and JUANES R., *Advances in Water Resources*, **144** (2020) 103684.

[21] THIELE U., *Colloids and Surfaces A*, **553** (2018) 487.

[22] BELTRAME P. and THIELE U., *SIAM J. Appl. Dyn. Syst.*, **9** (2010) 484.

[23] BELTRAME P., KNOBLOCH E., HÄNGGI P. and THIELE U., *Phys. Rev. E*, **83** (2011) 016305.

[24] MITLIN V., *J. Colloid Interface Sci.*, (1993).

[25] THIELE U., BRUSCH L., BESTEHORN M. and BÄR M., *Eur. Phys. J. E*, **11** (2003) 255.

[26] DICARLO D. A., *Water Resources Research*, **40** (2004) W04215.

[27] BAUTERS T. W. J., STEENHUIS T. S., PARLANGE J.-Y. and DICARLO D. A., *Soil Science Society of America Journal*, **62** (1998) 1185.

[28] FETTER C. W., *Applied Hydrogeology: Fourth Edition* (Waveland Press) 2018.
## 2.6 Conclusion

The water transfer model proposed in this chapter reproduces several water transfers in the presence of hydrophobic matter: water trapping, fingering and Darcy flow. Taking into account the free energy of the soil/water system improves our understanding of water transfer in porous media. In addition, taking repellent energy into account is fundamental for introducing the hysteresis phenomenon. During water transfer, there are multiple states of equilibrium varying according to the average saturation between partial wetting and total wetting. In other words, the dynamics change according to the average saturation of the soil, particularly during the wetting/drying cycle.

This model offers a new way of representing water transfer in porous media in presence of hydrophobic matter. Its general expression can be adapted to different porous media, such as soil in the presence of amphiphilic matter. This is a key point in the modeling, as in the soil physics. However, the promised results must be compared with experimental data, particularly in the presence of amphiphilic matter. It is important to define the model of the porous medium studied, i.e. to define the physical parameters of the model. In other words, we need to characterize the porous medium in the presence of an amphiphilic substance. This study is presented in the context of the rhizosphere in the part III.

In the next chapter, we propose to present the numerical tool needed to simulate the water transfer present in this manuscript.
# Chapter 3

## Development and use of numerical tools to study water transfers

The aim of this chapter is to present the numerical models developed and used to study our fourth-order equation introduced in chapter 2. These modeling tools are used to define temporal solutions and to study dynamic behavior as a function of model parameters.

These tools are fundamental in order to numerically analyze the behavior of our fourth-order degenerate parabolic equation. Since analytical study is extremely difficult due to the strong non-linearity present in the equation, numerical methods are required to study water transfer in soil in the presence of amphiphilic matter.

In this chapter, we introduce the basics of numerical modeling: to present the modeling difficulties for the same type of equation in the literature, and to justify the development of a model with the finite element method during this thesis (section 3.1). To facilitate the reader's understanding of the developments, we present the principle of the finite element method in section 3.2 and the concept of spatial and temporal meshes in section 3.3. Then, we present and discuss algorithms developed to study the temporal dynamics of water transfer in a porous medium in section 3.4. Finally, we present bifurcation tools using continuation algorithm with the mean saturation as bifurcation parameter in section 3.5.

### Contents

<table><tr><td><strong>3.1</strong></td><td><strong>State of the art of the time integration</strong></td><td><strong>78</strong></td></tr><tr><td>3.1.1</td><td>Numerical methods</td><td>78</td></tr><tr><td>3.1.2</td><td>Numerical time schemes</td><td>79</td></tr><tr><td>3.1.3</td><td>Numerical problem statement</td><td>79</td></tr><tr><td>3.1.4</td><td>My approach</td><td>81</td></tr><tr><td><strong>3.2</strong></td><td><strong>Finite element method</strong></td><td><strong>82</strong></td></tr><tr><td>3.2.1</td><td>Weak formulation</td><td>82</td></tr><tr><td>3.2.2</td><td>Galerkin method</td><td>82</td></tr><tr><td><strong>3.3</strong></td><td><strong>Meshing</strong></td><td><strong>84</strong></td></tr><tr><td>3.3.1</td><td>Spatial mesh</td><td>84</td></tr><tr><td>3.3.2</td><td>Time mesh</td><td>87</td></tr></table>
<table><tr><td><strong>3.4 Numerical temporal scheme</strong></td><td><strong>91</strong></td></tr><tr><td>3.4.1 Implicit Crank-Nicolson scheme</td><td>92</td></tr><tr><td>3.4.2 Implicit Euler scheme</td><td>92</td></tr><tr><td>3.4.3 BDF2 scheme</td><td>94</td></tr><tr><td><strong>3.5 Numerical Bifurcation Analysis</strong></td><td><strong>95</strong></td></tr><tr><td>3.5.1 Principle of continuation method and bifurcation analysis</td><td>95</td></tr><tr><td>3.5.2 Implementation in PDE2PATH</td><td>97</td></tr><tr><td><strong>3.6 Conclusion</strong></td><td><strong>100</strong></td></tr></table>

# 3.1 State of the art of the time integration

## 3.1.1 Numerical methods

Numerical methods are chosen and analyzed according to different mathematical properties: convergence, stability, consistency, etc. They enable the exact solution of continuous equations to be linked to the exact solution of discretized equations and to the numerical solution obtained. In addition to these mathematical criteria, there are computing criteria: CPU costs, accuracy, simulation time and memory used.

In the field of applied mathematics, many mathematical tools are used to numerically approximate the solution of complex problems (wave dynamics, mechanics, thermodynamics, ...) modeled by partial differential equations (PDE). There are three main numerical methods for describing the solution of a partial differential equation from a continuous problem into an approximate discrete problem: Finite Difference [234], Finite Element [123, 225, 81, 46] and Finite Volume [99].

Finite difference methods are based on the development of Taylor series (derivative approximations). They are simple to construct, and form the basis of most industrial and research codes in computational fluid mechanics. Their limitation to simple geometries and the consideration of boundary conditions, in particular Neumann boundary conditions, are their main disadvantages. As a result, these methods are mainly used to study the time evolution of the numerical solution.

As opposed to Finite Difference methods, Finite Volume and Finite Element methods are based on approximations of integrals. Finite Volume methods are based directly on the strong form of the equation to be solved. They evaluate the exact expressions of the mean value of the solution on a given volume. These data is used to construct approximations of the solution.

Finite Element methods are based on a variational formulation of the equation which is considered as the weak form. These local approximations of the solution, using local data, form a global approximation by joining them together. Both methods can be easily applied to unstructured meshes, and are generally better suited to the numerical sapce treatment of the problems.

Finite Volume methods have been developed for conservation laws and do not take into account the complexity of the mesh in discretizing conservation laws. However, the geometric characteristics of the mesh can influence the dynamics, for example when diffusive flows are considered. Conversely, Finite Element methods depend on the geometry of the domain and using mixed formulations based on mathematical theorems [47], they allow to take into account local and global physical dynamics [123].
### 3.1.2 Numerical time schemes

Numerical methods for solving problems are based on numerical schemes. There are two types of numerical scheme:

* Explicit schemes based on solutions at present and past times, i.e.

$$ \dot{u} = f(u_n, u_{n-p}), \qquad (3.1) $$

with $u$ the unknown, $f$ a continuous function and $u_n$ the solutions at each instant $t_n$. Note that $t_n$ indicates the present time and $t_{n-p}$, with $n \ge p > 0$, are the past instants.

* Implicit schemes that are dependent on the future moment, i.e.

$$ \dot{u} = f(u_{n+1}, u_{n-p}), \qquad (3.2) $$

with $u_{n+1}$ the solution at each future instant and $n \ge p \ge 0$.

Implicit schemes are unconditionally stable in the linear case, i.e. they do not develop numerical instability except in the case of non-linearity. However, explicit schemes are numerically unstable if we do not impose a strong condition on the time step $\Delta t$ (time duration) similar to the Neumann's stability criterion [203].

Examples of single-step and multi-step numerical schemes are given in Appendix A.

### 3.1.3 Numerical problem statement

Eq. (2.40) defined in chapter 2, can be written,

$$ \phi \partial_t s = -\nabla \cdot [K(s) \nabla (\Pi_w(s) + \Gamma \Delta s + \rho g z)] = F(s, x). \qquad (3.3) $$

This fourth-order partial differential equation is a degenerate parabolic equation [158] representing water transfer in a porous medium. It is similar to the lubrication [257, 236], Cahn-Hilliard [349, 152] or Kuramoto-Sivashinsky [193] equations. The strong non-linearity of these equations makes it impossible to define an analytical solution. Therefore, we need to develop numerical methods to solve this problem numerically.

The Navier-Stokes equations, classically used in fluid mechanics, are solved numerically using numerical methods. However, these equations are of second order, whereas our equation is strongly non-linear and of fourth order, i.e it is more difficult to design numerical methods to simulate the physical behavior of the liquid. Tuckerman et al [339] has shown that numerical instabilities arise from higher-degree differential operators. For example, for Navier-Stokes equations, it's the Laplacian, and for lubrication equations, it is the Bilaplacian [22]. Both are diffusion-related operators.

Tuckerman et al [339] proposed to work with a semi-implicit scheme separating the linear operator $\mathbf{L}$ from the nonlinear $\mathbf{N}$ part, i.e. $F(s, x) = \mathbf{L}s + \mathbf{N}(s, x)$. This scheme is stable if $\mathbf{L}$ is the Jacobian matrix $\mathbf{J} = \mathbf{J}_{s_0} = D_s F(s_0, x)$. In order to speed up linear system resolution in the semi-implicit method, it is convenient to take for linear part the higher degree differential operator instead of $\mathbf{J}$.
This operator contains very large negative values as the the eigenvalues of **J**. Treating this operator implicitly will ensure the stability. Applied to the case of the Navier-Stokes equation, the linear operator **L** is the Laplacian operator, and for the lubrication and for the Cahn-Hilliard equations, **L** is the Bilaplacian operator.

Nevertheless, Beltrame and Thiele [22] have shown that this method does not succeed for the lubrication equation whenever wetting is partial, due to large variations in the disjoining pressure in the neighbourhood of the wetting front. It is important to treat implicitly the whole linear part, i.e. **J**, and not only the Bilaplacian. However, solving linear system involving the JAcobina matrix raise new problems. Moreovre there is no general efficient method to perform the inversion. They propose to use an exponential propagation time integration scheme [137, 175, 331] combined with a Cayley-Krylov projection [359] and an incomplete LU decomposition. This method has proved stable and allows the application of an adaptive time step [36, 22]. However, this numerical development, implemented in finite difference, depends on an arbitrary parameter during the Cayley-Krylov projection, which influences the simulation time. If this parameter is not adapted to the problem, it can also lead to instability due to spurious modes. Finally, Finite Difference method is not suitable to adaptive mesh. At the wetting front, a variable mesh size is crucial to avoid numerical instabilities. A Finite Element approach would be better for the spatial discretization.

In the literature, other Finite Difference numerical schemes have been proposed to integrate this type of equation: Chebychev pseudo-spectral method [85] or Positivity-Preserving Numerical Schemes [373]. However, the Chebychev pseudo-spectral method proposed by Cueto and Juanes [85] does not take into account a non-monotonic disjoining function at the origin of numerical instability, and the Positivity-Preserving Numerical Schemes proposed by Zhornitskaya and Bertozzi [373] does not take into account a disjoining term. As a result, these models do not suffer from the problems mentioned above for the development proposed by Beltrame and Thiele.

The study of numerical schemes without the disjoining term allowed to develop a new numerical method for solving Eq. (3.3). At the same time as Zhornitskaya and Bertozzi, Grün and Rumpf [160] proposed a similar version of Positivity-Preserving Numerical Schemes by considering a Finite Element scheme. Their scheme is based on entropy estimation (in one space dimension [26, 29], in multiple dimensions [263]). They based their scheme on conservation of the positivity of the discrete solution using an entropy formulation. When $\Pi_w = 0$ and $g_r = 0$, Eq. (3.3) admits the existence and non-negativity of solution (one space dimension [28] and higher space dimensions [158, 120]). These models are considered as precursor notions for stronger regularity results as well as results on finite propagation velocity (in one space dimension [27], in several space dimensions [30]) or on a waiting time phenomenon [264].

Historically, Barrett et al are pioneers of a finite element scheme guaranteeing non-negativity of solutions (for the case $\Pi_w = 0$ [11] and in the context of the Cahn-Hilliard equation [12]). At each time step, they solve a second-order elliptic variational inequality by imposing non-negative solutions. However, this solution is not necessarily as regular as the limit of discrete solutions obtained in [160], as the discrete Laplacian operator of the discrete solution is not controlled.

Later, Grün and Rumpf [161] adapted this numerical scheme to Eq. (3.3). The addition of this strong non-linear term introduces numerical instabilities. Grün and Rumpf propose to solve this
problem by decomposing the disjoining function $\Pi_w$ such that

$$ \Pi_w = \Pi_w^+ + \Pi_w^-, \qquad (3.4) $$

where $\Pi_w^-$ is the destabilizing function and $\Pi_w^+$ is the stabilizing function. They propose to explicitly treat only the destabilizing function $\Pi_w^-$. This strategy is similar to the method used for discretizing nonlinear source terms in semi-linear second-order parabolic equations [121]. Note that the primitive of the destabilizing function is concave, while the primitive of the stabilizing function is convex and non-negative. Naturally, this remark applies to the free-energy density presented in chapter 2. Thus, their numerical scheme achieves a strong convergence [159] in 1D of the discrete solution and simulation of dewetting in agreement with physical experiments [161]. Note that the method developed by Grün and Rumpf is efficient, but requires a good mesh to limit the degeneracy of the solution due to mobility.

### 3.1.4 My approach

We opt for the development on a Finite Element method to numerically solve the equations (2.40a) and (2.40b). This choice is based on the fact that this method is supported by mathematical theories such as the concepts and properties of functional analysis, such as Hilbert and Sobolev spaces [47, 46] and the spatial adaptive meshing criteria. Moreover, the non-linearity of the problem makes it very difficult to find a strong solution directly, then that justifies to not use the method of Finite Volumes and Finite Differences. Instead, we will use a time scheme based on Finite Difference methods such as multi-step methods (Appendix A).

Note that we need to develop this numerical model in a high-performance software program that allows to use finite element methods. We have therefore studied several software options: FreeFem++, Matlab, COMSOL or Vsoil.

One of the main issues is the software's handling of non-linearities. For this reason, we quickly excluded COMSOL and FreeFem++. The commercial software COMSOL does not handle the non-linear form of our equations, and it is very expensive. In the case of FreeFem++, the advantage is that it is based on finite elements and it is free, but the implementation of the non-linear parts was not intuitive.

Then we considered an implementation using the Vsoil software developed at UMR EMMAH, but the development was limited to 1D simulations and did not allow the development of adaptive meshes. Therefore, it was decided to use Matlab. This payware was used to develop Beltrame and Thiele [22] Finite Difference method. Moreover, code for post-processing of data is easily implemented for analysis of results with this software. It should be noted that the FEM Operations Toolbox is available for a fee. However, the tools available in the toolbox are not always suited to our non-linear elementary matrices. Additional numerical methods had to be developed for the numerical resolution of Eq. (2.40a) and (2.40b).

We will return later in the section 3.4 to the development of a more powerful code in particular developed by Hartmann and Thiele [323, 165] on similar equations. This development is implemented in C++ using the oomph-lib [169] library developed and maintained by Matthias Heil and Andrew Hazel (School of Mathematics, Oxford Road, University of Manchester).

We propose to review some of the basics of finite element methods (section 3.2) and spatial and
temporal mesh construction (section 3.3), in order to facilitate the reader's understanding of the construction of the numerical methods developed to solve the equations (2.40a) and (2.40b) presented in section 3.4.

## 3.2 Finite element method

The finite element method is a numerical method for solving partial differential equations with a weak or variational formulation in a discrete domain $\Omega_k \subset \Omega$, where $\Omega$ is a domain of $\mathbb{R}^N$ representing the soil studied. An algorithm is implemented to define an approximate solution of this PDE on a compact domain with conditions on the edges and/or in the compact interior.

### 3.2.1 Weak formulation

A statement of a physical problem governed by PDE, generally expressed in strong form, is written in weak form, called weak formulation or variational formulation. Let define the mathematical space needed to express this formulation.

Let $\Omega$ the soil studied without boundaries and $\bar{\Omega}$ the whole soil studied during a period of time $[0, T]$, with $T \in \mathbb{R}$. Let the couple $(s, \frac{\delta F}{\delta s}) \in C^2(\bar{\Omega} \times [0, T]) \times C^2(\bar{\Omega} \times [0, T])$ a classical solution of Eq. (2.58) and $V_1$ and $V_2$ two subspaces of $H^1(\Omega \times (0, T))$. Note that, subspaces $V_1$ and $V_2$ depend on the boundary conditions (see section 2.2.3).

Eq. (2.58) is related to the following weak formulation:

$$ egin{cases} 
\text{Find } (s, \frac{\delta F}{\delta s}) \in V_1 \times V_2 \text{ such that,} \\ 
\int_{\Omega} \partial_t sv = \int_{\partial\Omega} K(s) \partial_n \frac{\delta F}{\delta s} v - \int_{\Omega} K(s) \nabla \frac{\delta F}{\delta s} \nabla v \\ 
\int_{\Omega} \frac{\delta F}{\delta s} w = \int_{\Omega} \partial_s f_w w - \int_{\partial\Omega} \Gamma \partial_n s w + \int_{\Omega} \Gamma \nabla s \nabla w - \int_{\Omega} \rho gzw, 
\end{cases} 
\quad (3.5) $$

where $v$ and $w$ are test function.

This formulation, based on properties of functional analysis, is used to approximate the solutions of the problem in weak form. To simulate this model, we need to determine a discrete problem.

### 3.2.2 Galerkin method

In numerical analysis, Galerkin methods are a class of methods for transforming a continuous into a discrete problem. The continuous problem is given by the abstract formulation:

$$ egin{cases} 
\text{Find } u \in V \text{ such that} \\ 
a(u, w) = f(w), \quad \forall w \in W, 
\end{cases} 
\quad (3.6) $$

where $V$ and $W$ are Sobolev spaces, $a$ is a bilinear form defined on $V \times W$ and $f$ a linear form on $W$. Note that $a$ and $f$ are approximated to transform them into linear form if they are non-linear.
Because the problem Eq. (3.5) is a non-linear system, its linear formulation is required in order to define the approximated solution $u_h$:

$$ \left\{ \begin{array}{l} \text{Find } u_h \in V_h \text{ such that} \\ a_h(u_h, w_h) = f_h(w_h), \quad \forall w_h \in W_h, \end{array} \right. \qquad (3.7) $$

where $V_h \subset V$ and $W_h \subset W$ are Sobolev subspaces, $a_h$ is a bilinear form defined on $V_h \times W_h$ and $f_h$ a linear form on $W_h$.

In our case, the abstract formulation is defined by linearising the non-linear functions $K$ and $\partial_s f_w$. Once the system is linearised as in section 2.3.3, Galerkin's method is used to approximate the solution $u$ by a solution $u_h$ using a discrete linear system. Therefore, the subspaces of our studied Sobolev spaces must be defined.

Let $N$ and $M$ denote respectively the dimensions of the subspaces $V_h$ and $W_h$, i.e. $N = \dim V_h$ and $M = \dim W_h$. In other words, the domain is discretised into a mesh of $N$ elements in space $V_h$ (section 3.3.1). The families of vectors $\{\varphi_1, \dots, \varphi_N\}$ and $\{\Phi_1, \dots, \Phi_M\}$ are the base elements of $V_h$ and $W_h$ respectively. The approximate solution $u_h$ is then decomposed onto each element from the basis $V_h$,

$$ u_h = \sum_{i=1}^{N} U_i \varphi_i. \qquad (3.8) $$

This decomposition associates to the discrete weak formulation a matrix system, derived from linear algebra, such that

$$ A(U^{-})U = b(U^{-}), \qquad (3.9) $$

where $U^{-}$ is a vector known from our linearized problem.

The elements of the matrix $A$ are given by $A_{ij} = a_h(\varphi_i, \Phi_j)$ with $i \in \{1, \dots, N\}$ and $j \in \{1, \dots, M\}$, the unknown is given by $U \in \mathbb{R}^N$ and the RHS by $b \in \mathbb{R}^M$ with $b_i = f_h(\Phi_i)$.

However, the matrix form obtained depends on the time scheme used to approximate the continuous system (Appendix B). In addition, the elements used to define these matrices influence the formulation of the matrix system. Here we choose to use Lagrangian elements to represent functions of the form $\varphi_i = \Phi_i$ for all $i \in 1, \dots, N$ with $M = N$. In other words, we use polynomial elements and resolve a matrix system with a square matrix $A$.

In our problem, Lagrange polynomials $P_1$ are at least used to have a correct approximation of the weak solution. These elements are linear function as displays in Figure 3.1. In our development, we propose to use Lagrange elements $P_2$ to obtain a higher accuracy of the weak solution without drastically increasing the computation time. The higher the degree of the Lagrange elements is chosen, the greater the construction time and complexity of the global matrices are. This criterion influences the time taken to solve the problem at each time iteration.

A weak solution is approximated numerically using the finite element method on each element of the discrete domain. A discrete solution needs to be defined using a method for discretizing the unknowns.
<figure><img src="image_50.png" /><figcaption>Figure 3.1: Quadrature function of test function $P_1$.</figcaption></figure>

## 3.3 Meshing

Meshing is an important element in the study of physical behavior, particularly when it is governed by a partial differential equation, as demonstrated in section 3.2 with the use of Galerkin method. It represents the decomposition of a continuous spatial or temporal domain into a discrete domain. In this section, we assume that spatial meshing 3.3.1 is associated with the finite element method and temporal meshing 3.3.2 with the finite difference method.

### 3.3.1 Spatial mesh

A spatial mesh is the discretization of a continuous space by finite and similar elements. The same category of geometric shape is used, but the proportions may be different. Elements are generally found in the form of a segment in 1D, or in the form of a quadrilateral or triangle in 2D. This discretization must respect certain mathematical rules to be well defined. Let a set $\Omega$ and $\{T_k\}$, with $k \in \{0, 1, \dots, m-1\}$, a family of geometric form of $m$ elements such that the discrete domain $\Omega_k$ is expressed in the form

$$ \Omega_k = \bigcup_{k=0}^{m-1} T_k. \qquad (3.10) $$

This discrete domain is a conformal mesh if two elements $T_k$ and $T_j$ with $k \neq j$ and $k, j \in \{0, 1, \dots, m-1\}$ satisfy one of this criteria:

* $T_k$ and $T_j$ are disjoint, i.e $T_k \cap T_j = \emptyset$;

* the intersection of these elements is a common edge or vertex.

Figure 3.2 shows some examples of meshes conformal and non-conformal for a rectangular set.
<figure><img src="image_51.png" /><figcaption>Figure 3.2: Meshing in a rectangular set</figcaption></figure>

These conformal meshes correspond to a tiling of the studied domain. The vertices of the elements of our tiling correspond to the nodes of the mesh. They define the discretization points that allow the continuous problem to be simplified into a discrete problem, as mentioned in section 3.2.2. Note that the approximate solution of the initial continuous problem is calculated at each mesh node. This discrete space can be associated, in our case study, with the experimental measurements obtained locally (in 1D) or in a detection zone (in higher dimension) by sensors in the soil. We can distinguish two types of node configuration, i.e. two categories of mesh:

* uniform meshes, consisting of a single geometric shape, i.e. element of the mesh has the same geometric shape;

* non-uniform meshes, made up of a family of geometric shapes in which at least one shape in the mesh has different proportions to another element.

## Uniform mesh

In the case of a uniform mesh, the link between the numerical mesh and the experiment is facilitated because we can place our sensors at constant intervals. However, in the study of the dynamics of water in a soil, this spacing can be a source of error in the numerical resolution. If the element size is too large, we risk to miss a precursor behaviour of the dynamics and to move the simulated dynamics far away from the experimentally observed dynamics. If the element size is too small, the computation time during increases significantly and due to the multiple approximations of the simulated dynamics errors are generated. Furthermore, from an experimental point of view, this represents a very high material cost, which is not always possible to set up in the field. Therefore, the uniform mesh is not the best choice for non-linear problem.
## Non-uniform mesh

The aim of non-uniform meshing is to strategically place the nodes of the mesh to better capture the behaviour of water in the soil. To achieve this, it is necessary to adapt their location using mathematical tools such as error estimator. These tools make it possible to estimate a priori the error made on each mesh element. This estimation attempts to adapt the error in the mesh using a posteriori estimation error [347, 350, 123]. An algorithm is implemented to optimise the a posteriori error on each element and improve the mesh. This algorithm consists of 4 main steps:

* step 1: construct an uniform mesh of $N = 8$ elements;

* step 2: estimate the a priori error on each element;

* step 3: adapt the mesh to optimise the error on each element;

* step 4: repeat steps 2 and 3 until to obtain a suitable mesh within the limit of $n_{iter} = 8$ remeshing.

We propose to use the Z2 error estimator in step 2 [378, 377]. This estimator evaluates the $\epsilon$ error committed with respect to the approximated solution on each element. The error is considered acceptable if it satisfies the following criterion:

$$10^{-5} \le \epsilon \le 10^{-3}. \qquad (3.11)$$

In the case where the error on an element does not satisfy the criterion, the mesh elements must be adapted in order to satisfy the criterion. Figure 3.3 displays the mesh adaption. If the error is too small, the element is unrefined to reduce the resolution time and accuracy (Figure 3.3a). If the error is large, then the element is refined to reduce the error (Figure 3.3b).

<figure><img src="image_52.png" /><figcaption>Figure 3.3: (a) Unrefined and (b) refined spatial mesh with triangle elements.</figcaption></figure>

By analogy with an experimental device, the new node position is not necessarily at the same position as the sensor initially placed on a sample. But it does give a better estimate of the numerical solution. Therefore, a non-uniform mesh provides a better reproduction of the dynamics, but is more difficult to compare with experience at a given position and time.

The mesh defined at a time $t_n$ is not necessarily adapted to the dynamics at the next time. It is important to be able to adapt the spatial mesh to the dynamics during the time evolution. However, the same approach as explained above should not be applied. Repeated adaptation of the domain discretization increases the calculation time and introduces numerical noise that can disturb the system. To limit these problems, we propose to generate an optimal mesh during the initialization of the problem by applying the algorithm described before, then to limit the adaptation to the maximum per time iteration, $n_{iter} = 1$.
Spatial meshes synthesis

Two types of spatial mesh can be distinguished: uniform meshes and non-uniform meshes. Uniform
meshes are easily developed and the comparison of numerical data with experimental data is natural.
However, they do not evolve over time, i.e. they depend on the choice of initial discretization, which
can affect simulated behavior when the phenomenon is complex. Non-uniform meshes, on the other
hand, reproduce more accurately the numerical behavior of the phenomena simulated. However, their
development and their comparison with experimental data is more difficult than uniform meshes.

Therefore, the spatial mesh is an important element in the estimation of the numerical solution,
but the type of mesh chosen depends on the constraints we require and the behavior studied.

3.3.2 Time mesh

A time mesh is a discretization of time by discretization points $t_n$, for all $n \in \mathbb{N}$. Each discretization point is separated by a time interval which we call the discretization step $dt_n$, where $dt_n$ is the space between the time at instant $t_n$ and the following instant $t_{n+1}$, i.e $dt_n = t_{n+1} - t_n$. As with spatial meshing, there are two categories of time meshing:

* uniform meshes where the discretization step is constant during the entire study, i.e. for all $n \in \mathbb{N}$,
  $$dt_n = dt$;

* non-uniform meshes, where at least one discretization step varies during the study, i.e. there is at least one step between time $t_{m+1}$ and $t_m$, such that $dt_m \neq dt_n$, for all $n, m \in \mathbb{N}$ and $n \neq m$.

Uniform mesh

Uniform meshes are easy to set up numerically. Once the discretization step is defined, the time
interval is constant. This time division ensures that exact time comparisons can be made between
numerical output and experimental data measured at constant time intervals. However, the choice of
discretization step must be selected strategically. In approximating the behaviour of the dynamics, the
discretization step can be a source of loss of information about the transfer in the soil. If *dt* is very
large, then there is a risk of changing the dynamics under study, as a fundamental behaviour of the
dynamics would be missed. If *dt* is too small, the number of iterations increases, which increases the
calculation time and the error of the approximate solution due to the time scheme. These increase can
affect dramatically the solution. Since the dynamics of water in the soil contain regimes on different
time scales, it is not appropriate to choose a uniform time mesh, but rather to adopt a mesh that
adapts to the dynamics, such as a non-uniform mesh.

Non-uniform mesh

Non-uniform mesh is associated to time adaptation which is based on the estimation and correction
of the time discretization step. Figure 3.4 displays the diagram associated to the time step resolution.
We consider a profile of the dynamics at a time $t_n$ and try to determine the profile at a time $t_{n+1}$ such
that $t_{n+1} = t_n + dt_n$. We assume that we have determined a time step $dt_n$ at the previous instant
$t_{n-1}$. We try to solve the problem on the discrete domain at time $t_{n+1}$, using an explicit numerical
scheme. We present here two methods to adapt the discretization step, described in Figure 3.4 by the
box *Estimation adapt*: constant prediction method and the prediction-correction method.

<figure><img src="image_53.png" /><figcaption>Figure 3.4: Diagram of the time step resolution</figcaption></figure>

The two methods have the same management when the resolution algorithm diverges. Since we
have not succeeded in defining a profile at time $t_{n+1}$, i.e. the step $dt_n$ given initially does not allow
us to converge, we need to correct it. We take the profile defined at time $t_n$ and repeat the resolution
to determine the profile at time $t_{n+1}$ by halving the time discretization step previously used, i.e.
$dt_n = \frac{dt_n}{2}$.

In the case where the resolution algorithm converges, we have succeeded in defining the profile of
the dynamics at time $t_{n+1}$. We try to predict the future discretization step $dt_{n+1}$ by means of an
estimate. The constant prediction method used in the Vsoil software [202] and other software suggest
considering that the future discretization step increase with respect to a constant. Generally, this
constant is selected at 10% as shown in Figure 3.5. This method adapts the time mesh but requires
multiple iterations to define a time instant if the constant is not correctly chosen.

Figure 3.6 displays the diagram of the discretization time step estimation using the prediction-correction method. The prediction-correction method is a numerical method for estimating the discretization step using mathematical criteria. Since the estimation of the time evolution of water in the
<figure><img src="image_54.png" /><figcaption>Figure 3.5: Diagram of the discretization time step estimation using the constant prediction method</figcaption></figure>

soil is governed by a numerical scheme, an error in the approximation of the solution is committed. This error is called the Local Truncation Error (LTE), denoted by $d_n$. To maintain the stability of the numerical scheme through time, we have chosen to use an implicit scheme. These schemes, known for their unconditional stability, commit an LTE of order $k$. In other words, the LTE is written as a function depending on the $k^{th}$ derivative of the variable $u$ at a time $t_n$, $u_n^{(k)}$. This error is expressed as a function of the difference between the approximate solution defined by the implicit scheme $u_{n+1}$ and an explicit scheme of the same order $u_{n+1}^P$. For example, an implicit Euler scheme (BDF1) with an explicit Euler scheme, a second order of the Backward differentiation formulas (BDF2 scheme) with a Leapfrog scheme or a Crank-Nicolson scheme/Trapezoid rule with a second order of the Adams-Bashforth methods (AB2 scheme). More details are given in Appendix A. This trick makes it possible to express the LTE without depending on the derivatives of $u$ and therefore without the need of a new numerical approximation.

The Root-Mean-Square Error (RMSE) norm of the LTE on each element of the spatial mesh is used with the convergence tolerance $\epsilon$ of the resolution algorithm to define the predicted time step $dt_c$ [317] (Appendix A),

$$dt_c = \left( \frac{\epsilon}{||d_n||} \right)^{\frac{1}{k+1}}.$$

As with the spatial mesh, we introduce a bounding criteria associated with our adaptation term estimate $dt_c$:

$$0.8 \le dt_c \le 4. \qquad (3.12)$$

If the adaptation estimate conforms to the criteria, then the future time step increases such that:

$$dt_{n+1} = dt_c dt_n. \qquad (3.13)$$

Otherwise, there are two cases:

* if the estimate is below the criteria, the rate of increase in the discretization step is set at the lower limit of the bounding set, i.e. $dt_{n+1} = 0.8dt_n$. This choice prevents the predicted time step from being reduced too much;
<figure><img src="image_55.png" /><figcaption>Figure 3.6: Diagram of the discretization time step estimation using the prediction-correction method</figcaption></figure>

* in the opposite case, the rate of increase in the discretization step is set at the maximum limit of the bounding set, i.e. $dt_{n+1} = 4dt_n$. This condition limits the growth of the predicted time step.

This prediction-correction method is very useful for optimising iterations and global calculation time by choosing strategic time discretization points.

## Time meshes synthesis

As for spatial meshes, there are two types of temporal mesh: uniform and non-uniform.

Uniform meshes are easily implemented in numerical methods. Note that experimental measurements of a physical phenomenon are generally performed at constant time intervals. Thus, the comparison between experimental and simulated data is well adapted to uniform meshes. However, this type of mesh is less suited to dynamics with time change.

However, the dynamic is better fitted with non-uniform meshes. They allow to define a faster resolution method for a final dynamic time, which reduces simulation time and the number of calculus iterations required. As a result, the error of the solution after several iterations will be lower than for a uniform mesh. However, the development of this method requires more implementation work and the selection of the method used to best fit the dynamics is needed (e.g. Figure 3.5 and Figure 3.6). It should also be noted that this type of mesh is less suitable for the comparison with experimental data than the uniform mesh.

Therefore, temporal meshes are extremely important to follow the dynamics of a numerical solution. In our situation, temporal meshes are fundamental to our numerical development to optimize the estimation of the simulated solution and better follow the time dynamics of the problem.
## 3.4 Numerical temporal scheme

During this thesis, three methods of resolution are developed to determine the time evolution of the saturation in a porous medium. These methods are based on adaptive time meshes according to the following numerical scheme:

i. implicit Crank-Nicolson scheme;

ii. implicit Euler scheme (BDF1) with an explicit Euler scheme;

iii. Backward differential formula (BDF2) scheme with a generalized Leapfrog scheme.

These methods are all included in the same global algorithm structure presented in Figure 3.7. It displays the diagram of the time evolution depending on the time step $dt$. The global algorithm

<figure><img src="image_56.png" /><figcaption>Figure 3.7: Diagram of the time loop.</figcaption></figure>

consists in determining the equilibrium state of the dynamics by studying the time step variations $dt_n$ between two time points $t_{n+1}$ and $t_n$, $dt = t_{n+1} - t_n$. A time loop is defined in order to determine the equilibrium state of water saturation in the soil. This loop is initialized for a time step $dt = dt_0$ defined at time $t_0$ and runs as long as the condition satisfies the following criterion:

$$10^{-12} \le dt \le 10^{12}. \qquad (3.14)$$

We have adopted the criterion, Eq. (3.14), to control the numerically solved dynamics and determine the equilibrium state. When the condition defined by Eq. (3.14) is not satisfied, we exit the iterative loop in one of two ways: the algorithm has reached the equilibrium state, or it has failed. We consider
that the algorithm failed, if $dt < 10^{-12}$. The dynamics is too fast and that the numerical scheme or the physical formulation is not suitable. In addition, it controls the decimals that the computer is able to interpret as non-zero. Generally, the value $10^{-12}$ is defined as a numerical zero by the computer. However, the algorithm reaches the equilibrium state of saturation in the soil, if $dt > 10^{12}$. The dynamics of the water transfer is very slow, so the dynamics is considered stationary. Note that the algorithm may continue to run and never reach the equilibrium profile if it is in a periodic regime. In this case, an *a posteriori* study is required, with a bifurcation study (section 3.5) for example.

When the condition defined by Eq. (3.14) is satisfied, the algorithm has not reached equilibrium and we continue the resolution described in the Time step resolution diagram shown in Figure 3.4 until the loop condition is no longer satisfied.

The three versions of this algorithm are carried out according to the numerical schemes used (version i, ii and iii), which influences the resolution performance to integrate the dynamics between the initial time $t_0$ and the resulting equilibrium state. In the following, we present and discuss the performance of each version.

### 3.4.1 Implicit Crank-Nicolson scheme

In a first step, we developed the global algorithm (Figure 3.7) in Matlab, based on the algorithm block from Figure 3.4. To simplify the development of this first version, we used a uniform spatial mesh composed of the Lagrangian finite element $P_1$ and an adaptive time mesh according to Constant prediction with $adapt = 1.1$ (Figure 3.5). We used the Crank-Nicolson scheme (scheme of order 2) with the method described by Grün and Rumpf [161] in the algorithm block associated to the Solver Step in the Time step resolution diagram (Figure 3.4). In this way, we have implicitly treated the stabilizing part and explicitly the destabilizing part. The matrix formulation of our equation is defined using the implicit scheme and it is solved using an iterative resolution method (Figure 3.8): the Newton/Picard method (described in the Appendix D).

Therefore, without the Matlab Toolbox, we have constructed elementary matrices suitable for our matrix formulation. Note that the Matlab language is known for its vector calculations, which is an advantage for solving the weak formulation in matrix form. The iterative Picard/Newton method defines if the calculated solution converges or not after a maximum number of iteration $k_{itermax}$ for a given explicit solution $u_{n+1}^P$. In this version, we assume that $k_{itermax} = 20$ and the initialization of the iterative method is performed with the explicit Euler scheme.

This version of the algorithm reproduces several water transfer dynamics, such as intermittence in a stratified medium in 1D with periodic conditions, or horizontal diffusion in 1D with homogeneous Neumann conditions. However, these time simulations are extremely long to reach the equilibrium dynamics. In the case of the 1D horizontal diffusion, the simulation was approximately one day, which is not sustainable for more complex domains such as 2D.

### 3.4.2 Implicit Euler scheme

In this second version, we propose to adopt the structure of the first version, changing the adaptive method and the temporal scheme used. Thus, we have an implementation in Matlab, without using toolboxes and with a uniform spatial mesh using Lagrangian finite elements $P_1$.
<figure><img src="image_57.png" /><figcaption>Figure 3.8: Diagram of the solver step.</figcaption></figure>

In this version, we opt for the temporal adaptive mesh based on the prediction correction method described in Figure 3.6. This mesh is more adapted to the simulated dynamics and limits the divergent loop turns resulting from the iterative resolution of the Picard/Newton algorithm in the Solver Step block (Figure 3.8). The aim of this approach is to drastically reduce simulation time compared with the first version.

We used the implicit Euler scheme, also call 1-step Backward differential formula (BDF1; scheme of order 1), with the method described by Grün and Rumpf [161] in the algorithm block associated to the Solver Step in the Time step resolution diagram (Figure 3.4). The explicit Euler scheme (scheme of order 1) is used to define the initialization of the iterative resolution method and to predict the future time step $dt$. The prediction correction method is used to compensate the losing of accuracy due to the schemes of order 1 compare to the scheme selected in this version.

This implementation obtained the same dynamics as the first version but with a shorter simulation time. The simulation was four times faster, dropping from one day to 6 h. However, this improvement still does not make the simulation time sustainable for more complex domains such as 2D. It is needed to retain the adaptive prediction correction method to define the temporal mesh and to improve the optimization of the functions developed, and to switch to an adaptive spatial mesh.
### 3.4.3 BDF2 scheme

In this third version, we aim to improve on the performance of the second version. This involves optimizing the implementation of numerical tools and the development of certain tools such as the creation of an adaptive spatial mesh. The Matlab language is sufficiently limited to allow all its improvements to be completed during this thesis. Moreover, the creation of matrices adapted to an adaptive spatial mesh becomes very quickly complex. Therefore, we considered using a different programming language and new software to acquire its improvements. Since one of the future objectives after this thesis work is to be able to integrate this numerical development into the Vsoil platform, in which models can be implemented in Fortran or C++, we decided to switch to development in C++.

The third version has been implemented in C++, using the Oomph-lib library. This library, developed in Manchester by Matthias Heil and Andrew Hazel [169], contains a number of numerical tools that have been presented in the sections 3.2-3.3. In addition, these tools are suitable for gradient-type equations [236, 271, 321, 165].

We opted an adaptive spatial mesh using the Z2 error estimator composed of the Lagrangian finite element $P_2$. We also used an adaptive time mesh based on the prediction correction method described in Figure 3.6 with the 2-step Backward differential formula (BDF2; implicit scheme of order 2). The solution is defined by the Newton algorithm with the initialized solution $u_{n+1}^P$ given by the generalised Leapfrog scheme (explicit scheme of order 2). Note that this scheme is used to predict optimally the time steps in the adaptive time mesh during the prediction-correction method.

However, these schemes of orders 2 are two-step multi-step schemes. A preliminary initialization step must be added in order to run the algorithm. We have added this step according to the structure described in Figure 3.9. This preliminary step is essential for initializing the time loop in the case

<figure><img src="image_58.png" /><figcaption>Figure 3.9: Diagram of the global code.</figcaption></figure>

where the time scheme chosen is a multiple-step scheme of order two or higher. We have to define the preliminary step. It is necessary to choose, in the preliminary step, time multi-step schemes having the same order as the schemes used in the time loop, i.e. schemes of order 2. Note that the implicit and explicit Euler schemes are time multi-step schemes of order 1, so they cannot be used. We choose the implicit Crank-Nicolson scheme, also known as the Trapezoid rule, during the solver step in the
preliminary step and the Adams-Bashforth scheme (AB2) as the prediction of the first step of each iteration in the preliminary step. These two schemes are of order 2. Using these two schemes, we determine the solution at time $t_1$, for a given initial time $t = 0$, and estimate the next discretization step $dt_1 = t_2 - t_1$ using the method described in section 3.3.2 for the box *Estimation adapt* (Figure 3.6).

This implementation produced the same dynamics as the previous two versions, but with a much shorter simulation time (of the order of one minute). A 2D development on several domains has been established with this version of the code. In addition, we are able to handle different boundary conditions: Dirichlet, homogeneous Neumann, flux or periodic. This version has been used to produce the Figures 2.14. Note that a possible improvement to this version is to use the Grün and Rumpf method to improve convergence between the weak and strong solutions.

We have developed a numerical model to simulate the dynamics of water transfer in a soil in the presence of hydrophobic matter until an equilibrium state is reached. This result is limited to the initial parameter used for the simulation. It is important to be able to study these equilibrium states as a function of system parameters and to distinguish between different possible scenarios. A bifurcation study is a powerful tool to describe this behavior.

## 3.5 Numerical Bifurcation Analysis

Time integration of the system of Partial Differential Equations (PDE) governing dynamics is dependent on initial conditions and the computation time required to capture phenomena over a sufficiently long period. As a result, it provides a very limited view of possible phenomena, without really being able to interpret them. A more precise description of the dynamics in phase space can be obtained by using the continuation method or path-following method for stationary or periodic states in parameter space [22, 102]. The principle is to consider different parameters, where stationary [periodic] solutions are known, then, to follow these solutions and their possible bifurcations by continuously varying the parameters. Thus, for each parameter variation, the phase portrait is known more precisely. In particular, unstable directions of equilibria provide information on the behavior of the dynamics (stable, bounded, etc...) for regardless the initial conditions. Knowledge of these bifurcation scenarios of equilibrium branches provides an explanation for the emergence of possibly complex dynamics.

With the model we have developed, several asymptotic dynamics (steady-state or periodic dynamics) may coexist. Moreover, we pointed out the existence of thresholds, such as that fingering flow or water retained in a hydrophilic layer (see previous chapter). These asymptotic dynamics, called attractors, if they are stable, can be captured by the continuation method. Threshold determination corresponds mathematically to the existence of bifurcations. This theory is generally takes place in the mathematical theory of dynamical systems [374]. Although it is possible to track branches of periodic solutions [341] we focus only on stationary states. In the next sections, we present the continuation method applied to stationary state and the associate numerical scheme.

### 3.5.1 Principle of continuation method and bifurcation analysis

Let us consider a stationary solution $s_0$ of the dimensionless PDE given by Eq. (3.3)

$$ \partial_t s = \nabla \cdot (K(s, x) \nabla [-\Pi(s, x) - T \nabla s - g_r z]) = G(s, x, z, p), \quad (3.15) $$
where $p$ is a parameter as the $c$ concentration of EPS or $s_m$ the mean water content. Since $s_0$ is stationary, then

$$G(s_0, x, z, p) = 0. \qquad (3.16)$$

According to the implicit function theorem, Eq. (3.16) defines an unique function $s_0(x, z, p)$ as long as the jacobian matrix $\frac{DG}{Ds}$ is invertible. The function $s_0(x, z, p)$ is continuous w.r.t the parameter $p$. If we associate to the function $s_0(x, z, p)$ a number corresponding to its norm $||s_0(x, z, p)||$ then by varying continuously the parameter $p$, we construct a continuous curve $(p, ||s_0||)$: it is the continuation method. When the operator is non-invertible there is a bifurcation. Universal properties of these bifurcations occur regardless of the physical context. Either a turning point appears: it is the saddle-node bifurcation or one (or many) branch(es) emerges(emerge)/vanishes(vanish). The main bifurcations are the transcritical, pitchfork and the Hopf bifurcations. The pitchfork and the Hopf bifurcations appear in the study of chapter 2. The different branches with the bifurcation points constitute the bifurcation diagram.

The continuation method consists in tracking the stationary solution $s_0(x, z, p)$ when varying the parameter $p$. The used method is the pseudo-arclength method based on a tangent predictor - secant corrector method. This method allows to follow steady-states regardless their stability and allow the turning points (saddle-node bifurcations) to be passed (Fig. 3.10). The tangent predictor - secant corrector method is detailed below.

Let us assume the $s_0(x, p_0)$ stationary state known and it is not a bifurcation point. We seek the solution $s(x, p_0 + \delta p)$. According to Eq. (3.16),

$$\frac{DG}{Ds}(s_0, x, z, p_0) \delta s + \frac{DG}{Dp}(s_0, x, z, p_0) \delta p = 0. \qquad (3.17)$$

Note that in the discretized space, $\mathbf{J}_0 = \frac{DG}{Ds}(s_0, x, p_0)$ is the Jacobian matrix and $\frac{DG}{Dp}(s_0, x, z, p_0)$ is a column vector. Fixing a small variation $\delta p_t$ of the parameter $p$ and if $\mathbf{J}_0$ is invertible (no bifurcation) then

$$\delta s_t = -\mathbf{J}^{-1} \left( \frac{DG}{Dp}(s_0, x, z, p_0) \right) \delta p_t. \qquad (3.18)$$

The tangent predictor at $p_t = p + \delta p_t$ is deduced by

$$s_t = s_0 + \delta s_t. \qquad (3.19)$$

The variation $\delta p_t$ is done by in order to stay close to the steady-state:

$$||G(s_t, x, z, p_t)|| < \epsilon_t ||s_0||, \qquad (3.20)$$

where $\epsilon_t$ is a small parameter typically taken in the range $10^{-3} < \epsilon_t < 10^{-1}$.

Next, the secant corrector is computed using Newton's method to solve Eq. (3.16) close to $(s_t, p_t)$. The secant is taken orthogonal to the tangent $(\delta s_t, \delta p_t)$ to be able to follow the branch even in the neighborhood of a turning point. For the Newton iteration step that starts from $(s_t, p_t)$ the condition
reads

$$ \mathbf{J}_k s_{k+1} = -G(s_k, x, z, p_k) - \frac{DG}{Dp}(s_k, x, p_k) \delta p_k \quad (3.21) $$

$$ \langle u_t, \delta u_k \rangle = 0 \quad (3.22) $$

$$ s_{k+1} = s_k + \delta s_k \quad (3.23) $$

$$ p_{k+1} = p_k + \delta p_k \quad (3.24) $$

where $\mathbf{J}_k = \frac{DG}{D\delta s}(s_k, x, p_k)$ and $\delta u_k, \delta p_k$ are the unknowns at step $k+1$. Eq. (3.21) and the orthogonality condition (3.22) written as matrix equation are

$$ \underbrace{\begin{bmatrix} \mathbf{J}_k & D_p G(u_k, x, z, p_k) \\ u_t & 0 \end{bmatrix}}_{=\mathbf{N}} \begin{bmatrix} \delta u_k \\ \delta p_k \end{bmatrix} = \begin{bmatrix} -G(s_k, x, z, p_k) \\ 0 \end{bmatrix} \quad (3.25) $$

The continuation step requires inversions of the Jacobian matrix $\mathbf{J}$ (tangent predictor) and of the matrix $\mathbf{N}$ (Newton corrector steps). As for the tangent predictor, linear systems are invertible except at bifurcation points.

<figure><img src="image_59.png" /><figcaption>Figure 3.10: Sketch of a continuation step using the tangent predictor-secant corrector method.</figcaption></figure>

### 3.5.2 Implementation in PDE2PATH

Inverting the matrix raises the same numerical difficulties as those encountered with the time integration scheme. The challenge lies in filtering out unwanted eigenvalues with a large negative real part due to the presence of differentiation operators. Using Finite Elements method, a continuation code has been implemented in Beltrame and Thiele [22]. The authors employed Krylov reduction using the Cayley transform and an incomplete LU factorization to perform the continuation step and the bifurcation analysis, i.e. the computation of leading eigenvalues with the larger real part. This code was used at the beginning of the thesis and notably to carry out the bifurcation diagrams in the article [20] presented in chapter 2.
The latter code is however not flexible regarding boundary conditions. Moreover, it does not benefit from adaptive meshing as is possible with finite elements. Therefore, we propose an alternative algorithm based on the same continuations methods, but using finite element methods to discretize the spatial domain as in the temporal case. This development is carried out in Matlab using the pseudo-arclength [197, 122] continuation method available in the PDE2PATH [342] package. Solution branches are tracked in parameter space via tangent predictor and Newton corrector while varying the control parameter.

To illustrate the implementation steps, we consider the case of stratified media treated in chapter 2, where $\Omega$ domain is periodic (periodic boundary conditions) and the bifurcation parameter is the mean saturation $s_m$. Let us introduce the variable $w$ such that $w = s - s_m$, Eq. (3.16) reads

$$0 = \nabla \cdot (K(w + s_m, x)\nabla[-\Pi(w + s_m, x) - T\nabla w]). \quad (3.26)$$

In order to facilitate the implementation of this algorithm, we use the same decomposition method of the problem than temporal study by introducing the variable $\eta$ (analogous to the functional derivative of the free energy in the time study), i.e. Eq. (3.26) is rewritten,

$$\begin{cases} 0 = \nabla \cdot (K(w + s_m)\nabla\eta) - \nabla \cdot (K(w + s_m)g_r e_z), \\ \eta = -\Pi(w + s_m) - T\nabla w. \end{cases} \quad (3.27)$$

To have the uniqueness of the solution to Eq. (3.27) there are usually additional conditions. In our case, we must impose the conservation of matter $\int_{\Omega} w = 0$. Moreover, in the presence of a continuous symmetry, i.e. in our case invariance along $z$, there is a continuum of solutions translated in the $z$-direction. To fix a solution $w(x, z)$, we need to prohibit variations $dw$ that induce translation wrt $z$, namely such that $w(x, z) + dw = w(x, z + dz)$. Since the translational symmetry group wrt $z$ is generated by the operator $\partial_z$, the solution must be orthogonal to $\partial_z w$.

In a similar manner to the Lagrange multipliers, we add to Eq. (3.27) two parameters $\lambda$, resp. $v$, associated to the mass conservation, resp. translation symmetry, such as

$$\begin{cases} 0 = \nabla \cdot (K(w + s_m)\nabla\eta) - \nabla \cdot (K(w + s_m)g_r e_z) - \lambda - v\partial_z w, \\ \eta = -\Pi(w + s_m) - T\nabla w. \end{cases} \quad (3.28)$$

These parameters act as continuation parameters. They are held at zero during all continuation steps, which means that the additional conditions are checked.

Using the finite element $\phi$ of Lagrange $P_1$ on the discretized weak formulation of Eq. (3.28), we obtain the matrix formulation $AU = b$ such that

$$\underbrace{\begin{bmatrix} 0 & \mathbf{K}_1 \\ -T\mathbf{K} & 0 \end{bmatrix}}_{\mathbf{A}} \underbrace{\begin{bmatrix} w \\ \eta \end{bmatrix}}_{U} = \underbrace{M_0 F(U) - v\mathbf{K}_z U - g_r \mathbf{K}_z K(w + s_m)}_{b}, \quad (3.29)$$
with the following complementary matrix forms:

$$F(U) = \begin{bmatrix} -\lambda & & \\ & -\Pi(w + s_m) - \eta & \\ & & \end{bmatrix}, \\ \mathbf{M}_0 = \begin{bmatrix} \mathbf{M} & 0 \\ 0 & \mathbf{M} \end{bmatrix}, \\ \mathbf{K}_1\eta = \int_{\Omega} K(w + s_m)\nabla\eta\nabla\phi.$$

Note that **M** and **K** are the mass and stiffness matrices respectively.
Finally, we use the L² norm to represent bifurcation diagrams

$$||s - s_m|| = ||w|| = \sqrt{\int_{\Omega} w^2(x, z)dxdz.} \qquad (3.30)$$

This norm reads in discrete form,

$$||s - s_m|| = \sqrt{\sum_{i} \sum_{j} \mathbf{M}_{i,j} w_i^2}. \qquad (3.31)$$

<figure><img src="image_60.png" /><figcaption>Figure 3.11: Horizontal saturation profile in a heterogeneous hydrophobic medium for different values of mean saturation $s_m$. Solid blue [dashed red] lines indicate stable [unstable] steady-state solutions. The horizontal line corresponds to solution branches with uniform distribution. Parameters $T = 0.001$ and $n = 2.68$.</figcaption></figure>

Figure 3.11 displays the bifurcation diagram of water distribution $s(x)$ steady-states in a homogeneous horizontal layer og length $L$ obtained using PDE2PATH. The soil model is given by Eq. (2.40) of the chapter 2 and the permeability an pressure functions are given by Eq. (2.59)-(2.62). Boundary conditions are periodic at $x = 0$ and $x = L$.

The trivial solution is the homogeneous distribution $s \equiv s_m$ displaying by the horizontal straight line
in Figure 3.11. we numerically calculate the linear stability of the corresponding these steady-states.
We detect two subcritical pitchfork bifurcations (see red dots on the horizontal line in bifurcation
diagram Figure 3.11). The bifurcated branch corresponds to an non-uniform water distribution. Near
the bifurcation points, the branch is unstable but becomes stable via a saddle-node bifurcation (red
dot in Figure 3.11). The solution profile on this stable branch corresponds to the coexistence of dry
and wet region as shown Fig. 2 of our publication [20]. The resulting bifurcation diagram is similar to
Fig. 3 in [20].

3.6 Conclusion

We have presented the numerical tools developed and implemented to study the temporal dynamics of water transfer and the equilibrium dynamics of solutions. These developments are based on the finite element spatial discretization method. These tools, used on the Eq. (3.3), allowed to highlight different types of flow in a homogeneous then stratified medium in the chapter 2: intermittence, fingering, water trapping.

The numerical tools developed for the temporal study are efficient for the case studies in this
thesis. We have the tools to manage a wide range of problems in several spatial dimensions (1D and
2D classical or axisymmetric) with various boundary conditions (Dirichlet, homogeneous Neumann,
periodic and flux). Note that this numerical development works for initial analytical profiles, but can
also be applied to experimental data from a file. This offers greater flexibility of execution. Moreover,
the numerical tools we have developed can be used to solve problems coupled with other dynamics,
as we will see in chapter 5. Additionally, we can improve this temporal development by using the
algorithm proposed by Grün.

Conversely, the continuations tools developed are limited to periodic cases. An improvement would
be to develop a version for different boundary conditions that could be applied to 1D and 2D do-
mains [138, 176, 139]. This version would help us to study a wider range of situations and to better
understand certain physical behaviors.

The development of these tools is important for studying the dynamics of water transfer in soils in
the presence of amphiphilic matter. In the following, we propose to focus on the application of these
tools to the rhizosphere.
# Part III

## Application to the rhizosphere
<figure><img src="image_61.png" /></figure>
# Chapter 4

## Characterization and parameterization of soil water properties in the presence of amphiphilic matter

The aim of this chapter is to characterize a soil properties in the presence of amphiphilic matter in order to parameterize the hydraulic properties required for water transfer modeling (part II). This characterization will be treated with reference to specific soil zone: the rhizosphere.

Soil wettability data in relation with amphiphilic matter are rather scarce in the literature. To calibrate our water transfer model applied to rhizosphere, we carried out complementary laboratory experiments in order to get more data. In this experiment, we will define and create an analogue rhizosphere soil, as an artificial porous medium, to get wettability tendencies and hydraulic behaviour in the soil.

In this chapter, we first present the design of the artificial rhizosphere medium with different mucilage concentrations: the selection of compounds for root mucilage and the porous substrate (section 4.1). The characterization of hydraulic properties is presented section 4.2, highlighting the wettability trends. Section 4.3 presents and discusses the parameters derived from the characterization of the model.

### Contents

4.1 An artificial rhizosphere medium 104

4.1.1 Choosing of root mucilage for building rhizosphere soil 104

4.1.1.1 Mucilage sampling from young germinated seed root: Aeroponic collection 105

4.1.1.2 Mucilage sampling from maize brace roots 106

4.1.1.3 Mucilage sampling from Chia seeds 107

4.1.1.4 Synthesis on root mucilage 109

4.1.2 The porous media design for rhizosphere analogue 109

4.1.3 Design and elaboration of the analogue rhizosphere soil 110

4.2 Characterization of hydraulic properties of the rhizosphere analogue medium 112
<table><tbody><tr><td>4.2.1</td><td>Materials and methods</td><td>113</td></tr><tr><td>4.2.1.1</td><td>Analogue rhizosphere sampling</td><td>113</td></tr><tr><td>4.2.1.2</td><td>Saturated soil hydraulic conductivity measurement</td><td>114</td></tr><tr><td>4.2.1.3</td><td>Soil water retention curve</td><td>115</td></tr><tr><td>4.2.1.4</td><td>Wettability and drop infiltration measurement</td><td>116</td></tr><tr><td>4.2.2</td><td>Results</td><td>118</td></tr><tr><td>4.2.2.1</td><td>Saturated hydraulic conductivity</td><td>118</td></tr><tr><td>4.2.2.2</td><td>Retention curve following single samples</td><td>119</td></tr><tr><td>4.2.2.3</td><td>Retention curve with the "one pressure - one sample" process</td><td>121</td></tr><tr><td>4.2.2.4</td><td>Wettability measurement and drop infiltration test</td><td>122</td></tr><tr><td>4.3</td><td>Parametrization of soil hydraulic properties accounting for variable wet-tability</td><td>124</td></tr><tr><td>4.3.1</td><td>Wetting function: disjoining pressure</td><td>124</td></tr><tr><td>4.3.2</td><td>Hydraulic conductivity</td><td>127</td></tr><tr><td>4.4</td><td>Conclusion</td><td>128</td></tr></tbody></table>

## 4.1 An artificial rhizosphere medium

The rhizosphere is the soil zone influenced by the root activity [171, 174, 259]. It is a biological hotspot which influences physico-chemical and hydraulic properties of this soil zone. In the rhizosphere, the roots release organic matter as rhizodeposits consisting of mucilage (ExoPolySaccharides - EPS, being a major component), border cells or simple molecules [185, 133], which highly stimulate microbial activity.

Microorganisms of the rhizosphere degrade the exuded root mucilage but also release their own Extracellular Polymeric Substances, generally different from those degraded [75, 255, 292, 4]. As we are looking for a porous medium with a mucilage concentration that limits microbial activity in order to maintain a stable mucilage concentration during the experiment, the artificial medium must be relatively simple to construct and low in nutrients.

This artificial medium must contain root mucilage that can be added with concentration to the rhizosphere.

### 4.1.1 Choosing of root mucilage for building rhizosphere soil

As stated above, in the soil, organic substances are exuded by the plant's roots, driving interactions with soil microflora and the rhizosphere effect [185]. A large quantity of these substances are produced in the form of hydrated mucilage. It is made up of between 15 and 90% of the plant's photosynthesis [71, 351, 185]. Measured into a sand medium, this secretion of mucilage into the rhizosphere is estimated between 0.05 and 50 mg of dry mucilage for a 1 g mass of dry sand or $5 \times 10^{-3}\%$ to 5% expressed in percentage of sand dry mass [376]. Furthermore, according to Landl et al [204, 206], the high concentrations of 1% would be found only near the tip and close to the root surface, while the lower concentrations would be located further behind the tip and away from the root surface.
Mucilage is a biohydrogel made up of EPS, proteins, minerals and lipids [247]. The proportion and chemical composition of these components vary according to the plant and its environment. However, it is estimated that the contribution of EPS is dominant. It amounts between 78% [246] and 94% [38] of the root mucilage.

EPS plays an important role in the biological and physic-chemical properties of the rhizosphere. It is carbon substrate and energy source for rhizosphere's microorganisms [195, 3, 243, 246] and influences water transfer in this zone [363, 89, 335, 65, 244].

Mucilage forms a gel by the cross-linkage of its polymer chains [283]. Once exuded, this gel interacts with the water and soil. The mucilage adsorbs to soil particles, creating a 3D mesh interconnecting the soil particles [38]. This mucilage link can facilitate cation exchange, water absorption and transfer in soil in the presence of mucilage [246]. When it comes to choose a representative molecule for mucilage analogue, embedding its various properties, synthetic molecules such as PolyGalacturonic Acid (PGA) were often used in 1990 and 2000 studies [148, 371, 13, 235, 266]. PGA was incorporated into synthetic porous media and calcium added for gel formation [376, 88, 13].

However, this substance has been debated in recent years, in particular because of its divergent physical effect compared with root mucilage [39]. In our experimental approach, we have therefore opted for the use of mucilage produced by plants. However, this needs to collect plant's mucilage and we tested different methods for mucilage collection.

### 4.1.1.1 Mucilage sampling from young germinated seed root: Aeroponic collection

Mucilage can be collected from young roots after germination. This sampling can be carried out using an aeroponics system and adapted for different plant species (Appendix E).

Aeroponics is a method of growing plants with roots hanging in air with a very high humidity level [376, 178].

We use the method described by Zickenrott [376] and Holz [376] (Figure 4.1a), with seedlings using a hermetically sealed box with two areas:

* in the bottom, a layer of demineralized water with aquarium bubblers;

* in the upper part, a metal grid tray where seeds are placed (about 100 – 500 seeds, depending on grain size).

After a few days in the system, the seeds germinate and grow seminal roots which start exuding mucilage near the tip, forming droplets with the high humidity level of the box.

During 2-3 days of growth, before roots reach the water layer, the mucilage is collected using a metal stick and a suction device. This suction device is made of a vacuum pump connected to a baby fly, also known as a nasal aspirator (Figure 4.1b).

The mucilage is accumulated at the extremity of the young root and collect each days. Collecting mucilage from the roots of seedlings on tray is rather quick, requiring less than half an hour per day. The collection is accumulated in the same vial during the 3 days of growth, for getting a sufficient quantity. With this low-cost method, the chemistry of mucilage collected from different plants can be studied (Appendix E). However, the collected amount (some mL) was not sufficient to produce our porous media impregnated with mucilage.
<figure><img src="image_62.png" /></figure>

<figure><img src="image_63.png" /><figcaption>Figure 4.1: Aeroponic mucilage collection system. (a) Aeroponic system. (b) Aeroponic collect of root mucilage.</figcaption></figure>

### 4.1.1.2 Mucilage sampling from maize brace roots

Another possibility to get plant root mucilage is to use roots that emerge in air, not in soil. Maize develops such roots, called brace roots, that emerge first in air and later go into soil. These brace roots at the base of maize stem exude mucilage at their tip [246]. This mucilage is similar to the one exuded by underground roots [71, 258].

Figure 4.2a shows maize brace roots secreting mucilage in the field. This photo was taken in a field at INRAE's experimental center in Avignon.

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

<figure><img src="image_65.png" /><figcaption>Figure 4.2: Collection of brace maize root mucilage (a) Brace maize roots in field with mucilage secretion (b) Brace roots maize mucilage collect in the lab.</figcaption></figure>

We adapted the protocol proposed in the literature [25, 2] to collect mucilage of maize brace root. We sampled maize plants in the field by cutting them at ground level in the field. These sampled plants were selected based on the following two criteria: high numbers of brace roots, and these roots do not touch the soil. After sampling, the plants were transferred to the laboratory. The roots were then cleaned to remove soil using demineralized water and a toothbrush. Cleaning the plant is important to minimize soil contamination of the mucilage that will be collected. After washing, the maize stem
base with roots is placed in a hydroponic solution with calcium-rich mineral water. Collection of mucilage takes place after 24 hours in the hydroponic system with the suction device during 2-3 days (Figure 4.2b).

Although this method gives access to root mucilage from field plants, (i) it necessitates a maize crop, at the good development stage, (ii) rather small amounts of mucilage are sampled: around 100 mL of fresh mucilage is collected over a period of several weeks using about a hundred maize plants (although this mucilage exhibits a high fresh/dry matter ratio: 0.93%). Moreover, the mucilage collected revealed traces of border cells and other soil deposits (Figure 4.3) and would necessitate further processing before use.

<figure><img src="image_66.png" /><figcaption>Figure 4.3: Brace root mucilage with soil and rhizodeposits.</figcaption></figure>

### 4.1.1.3 Mucilage sampling from Chia seeds

Chia seeds are seeds of *Salvia hispanica*, a flowering plant (from *Lamiaceae* family). Chia seeds produce mucilage on contact with water. Its mucilage has a chemical composition (primarily xylose, glucose, and uronic acids; [216, 3]) and physical properties such as viscosity similar to those of mucilage exuded by the roots of lupines (*Lupinus albus L.*) and maize (*Zea mays L.*) [282, 67, 3] and microbial Extracellular Polymeric Substances in soils [316, 59]. Its mucilage forms a gel that absorbs water when it is hydrated and becomes hydrophobic as it dries [3].

Because of the experimental problems of obtaining sufficient quantities of root mucilage [376, 198, 24, 39], chia seed mucilage has been chosen as a model mucilage since few years. It is easily and fastly available in large quantities. Collecting mucilage requires two steps:

1. mucilage production;

2. mucilage collection with a suction system (Figure 4.4).

For mucilage production, we placed 5 g of chia seed in a tube with 50 g of water and a magnetic stirrer shaked the mixture during 2 h [3]. During shaking, the seed's mucilage is exuded, gradually transforming the water in the tube into a gel. At the end of shaking, the fresh mucilage was ready for collection.

Collection is carried out using a suction system. The suction device is based on two soil filtration systems and a suction pump. Figure 4.4a is an illustration of the system developed (Figure 4.4b). Two
<figure><img src="image_67.png" /><figcaption>Figure 4.4: Chia seed mucilage collection system. (a) Sketch of the suction system. The red arrows represent the suction exercised by the pump located at the far right. (b) Photo of the collection system. (c) Chia seed mucilage collected after passing through the collection system. (d) Mixture of mucilage and chia seeds before (left) and after (right) passage through the collection system.</figcaption></figure>

soil filtration systems (A and B) are connected by a pipe. System B is connected by another pipe to a suction pump. Suction is imposed by the pump in both soil filtration systems and is represented by the red arrows in Figure 4.4a. Collection is carried out in device A, comprising two compartments separated by a filter with a mesh size of 750 µm. This filter separates the mixture of mucilage and seed deposited in the upper compartment. The mucilage is filtered and recovered in the lower compartment under pressure in the system (Figure 4.4c). A spatula is required to facilitate the filtration process. Chia seeds clog part of the mesh if not moved regularly. System B is used as a buffer to protect the pump from mucilage. A piece of fabric is placed in the lower compartment of this closed soil filtration device.

This collection method is very efficient and guarantees a huge quantity of mucilage: 20 mL of fresh mucilage is collected in 3 h time (2 h of mixing and maximum 1 h of aspiration). We note that we could collect almost all the fresh mucilage initially (Figure 4.4d). Only the mucilage near the seeds could not be collected. Indeed, the capillary forces binding the remaining mucilage to the chia seeds were too strong for our suction pump to extract it. We also observed that the fresh/dry matter ratio was correct (ranging from 0.4 to 0.6% according to our measurements and the literature [199, 178, 24]) for processing into artificial porous medium.

Therefore, chia seed mucilage is a candidate to be used as model mucilage for the creation of the rhizosphere analogue medium experiments.
#### 4.1.1.4 Synthesis on root mucilage

There are several methods for collecting root mucilage from plants at different stages. After testing some different methods, the mucilage amount obtained from most of these methods is too small for a further use in a rhizosphere porous medium analogue [376, 245]. Galloway [141] explain that was difficult to collect enough pure mucilage, as we have shown. Thus, collecting enough mucilage is a limiting factor for carrying out experiments on the hydraulic properties of this porous medium impregnated with mucilage especially when a number of replicates is needed. Within the panel of methods tested, the only acceptable candidate for mucilae production is from chia seed. This mucilage is similar to that of lupins and maize, despite the fact that seed- and root-derived mucilage differ in some aspects. However, this mucilage has specific characteristics that differ according to the stage of the plant and its environmental context [101].

In the following of this manuscript, we will consider chia seed mucilage as our model mucilage.

In the future, it would be relevant to carry out further studies, comparing the chemical and physical properties of this mucilage with other mucilages at different plant stages. An initial study of mucilage characteristics for different plant seedlings is presented in Appendix E. A comparison of these results with chia seed mucilage could be envisaged in order to estimate water behavior with another model mucilage.

### 4.1.2 The porous media design for rhizosphere analogue

Soils are extremely diverse, depending on pedogenetic processes at work, climate and vegetation/-land use and will present large range of properties (chemical, physical, but also biological), texture, structure,... This diversity will also impact rhizosphere properties, which itself is submitted to intense biological, chemical and physical interactions, and make quite complex the study of this specific medium. To examine the hydraulic properties linked to mucilage effects, we will thus opt for a "simple" medium in order to get a stable medium over time and to limit microbial action. Sand is often used as "simple" porous medium and is a good candidate for hydraulic experiments on an artificial rhizosphere.

<figure><img src="image_68.png" /><figcaption>Figure 4.5: Particle size of Fontainebleau sand from [112].</figcaption></figure>

Its structure remains unchanged in the presence of water, and as oligotrophic substrate, should limit microbial development. We used the "Fontainebleau" sand in our experiments, often used in
geophysics, geotechnical or porous media research. This sand is a pure silica sand (>99.5% SiO₂) made up of fine, round to sub-angular white grains [95, 34, 112, 6], with an uniform particle size distribution (mode and mean particle size of about 200 µm ± 100 µm [34, 109, 80] ; Figure 4.5). This substrate is classified as a very fine [215], which has an average dry density ranging from 1.39 to 1.77 g.cm⁻³ [34, 95]. This sand is available from suppliers (WRR). In order to control our medium, we selected a size range of particles between 200 and 250 µm by dry sieving using stainless steel sieves. Therefore, the homogeneous Fontainebleau sand we used has a particle size of between 200 and 250 µm. A parameterization of the hydraulic properties of Fontainebleau sand from different authors

<table><caption>Table 4.1: Hydraulic parameters of Fontainebleau sand. Parameters φ, θs, θr, Ksat, λ, Hcap are respectively the porosity, the saturated volumetric water content, the residual volumetric water content, the saturated conductivity, the shape parameter of the Brooks and Corey model and the capillary height. The symbol "X" corresponds that the data is not mentioned.</caption><thead><tr><th>Authors</th><th>φ</th><th>θ<sub>s</sub><br>[m<sup>3</sup>.m<sup>-3</sup>]</th><th>θ<sub>r</sub><br>[m<sup>3</sup>.m<sup>-3</sup>]</th><th>K<sub>sat</sub><br>[m.s<sup>-1</sup>]</th><th>λ</th><th>H<sub>cap</sub><br>[m]</th></tr></thead><tbody><tr><td>Doussan et al. [112]</td><td>0.36</td><td>0.36</td><td>0.013</td><td>3.09 × 10<sup>-4</sup></td><td>X</td><td>X</td></tr><tr><td>Bordes et al. [34]</td><td>X</td><td>X</td><td>X</td><td>5.69 × 10<sup>-5</sup></td><td>X</td><td>X</td></tr><tr><td rowspan="2">Allegre et al. [6]</td><td>0.36 – 0.37</td><td>0.355</td><td>0.11</td><td>1.72 × 10<sup>-4</sup></td><td>3.65</td><td>0.4</td></tr><tr><td>0.36 – 0.37</td><td>0.358</td><td>0.12</td><td>1.72 × 10<sup>-4</sup></td><td>3.88</td><td>0.45</td></tr></tbody></table>

is presented in Table 4.1. Note that the parameters are similar between the different studies. These values will be useful to compare with our control sand for retention measurements.

### 4.1.3 Design and elaboration of the analogue rhizosphere soil

#### General elaboration of the analogue rhizosphere soil

We chose Fontainebleau sand (section 4.1.2) and chia seed mucilage (section 4.1.1) as the main constituents of the analogue rhizosphere soil. The elaboration of these analogue media is carried out in several steps:

1. Mucilage collection;

2. Sand purification;

3. Wet/dry weight ratio of collected mucilage;

4. Mixing of sand / mucilage;

5. Sampling and humectation desiccation cycles.

The different steps of this design method are described in detail in Appendix F. Step 1 is the collection of a sufficient fresh amount of mucilage needed to build the medium, using the system described in section 4.1.1. In step 2 the sand is washed limiting the presence of metals and bacterial colonies. After completing these steps, the mucilage is mixed with the substrate. In this step, the quantity of mucilage to add per gram of dry soil is needed.

We designed our rhizosphere analogue medium with a range of concentrations that have been estimated in the rhizosphere. According to section 4.1.1, we decided to design media with concentrations
of 0, 0.1, 0.5 and 1 g of dry mucilage for 100 g of dry sand and expressed as percentage in the fol-
lowing (c = 1% correspond to 1 g of dry mucilage for 100 g of dry sand). Analogue media with higher
concentrations than 1% were not produced because of high amount of mucilage needed.

The wet/dry weight ratio of collected fresh mucilage (for each mucilage collection) was measured by oven drying (drying at 105°C for 24 h) and this ratio used in step 4 for mixing sand and fresh mucilage at the target concentration. In the step 5, samples are prepared for hydraulic experiments: saturated hydraulic conductivity, water retention and wettability measurement. During this step, humectation/desiccation cycles are performed to stabilize the structure of the medium. The humecta-tion/dessication method consists in humectating the medium by capillary rise from the bottom to limit the presence of air in the sample, then after saturation drying it at 30°C during 48 h [88, 199, 108].

Design of analogue rhizosphere soil

Two types of analogue rhizosphere soil were designed for the experimental samples. Samples format
are described in more detail in section 4.2.1.

Initially, undisturbed media were designed for our experiments. Figure 4.6 shows the different steps
required to elaborate analogue rhizosphere media. However, during drying of the sample in supports,

<figure><img src="image_69.png" /><figcaption>Figure 4.6: Preparation of undisturbed analogue rhizosphere soil.</figcaption></figure>

the surface of the medium was deformed: the medium rising to the edges of the sample support (Figure 4.7). In addition, we suspect a heterogeneity of the mucilage concentration in the medium: the surface in contact with the air seemed to have a stronger bond between sand particles than the lower part of the sample.

<figure><img src="image_70.png" /><figcaption>Figure 4.7: Photo of undisturbed analogue rhizosphere soil with c = 1%.</figcaption></figure>

We designed disturbed media to homogenize the mucilage concentration in all the samples. During the step 5, a wetting/drying cycle were realized on the sample before being mixed using a mortar and pestle (Figure 4.8). The mucilage adsorbed to the sand particles was uniformly distributed in the sample. During the new humectation phase, mineral water rich in calcium was used to promote the formation of bonds between sand particles by the mucilage. The presence of uronic acid in the mucilage
<figure><img src="image_71.png" /><figcaption>Figure 4.8: Preparation of disturbed analogue rhizosphere soil.</figcaption></figure>

in contact with calcium allows the formation of bonds with soil particles and ions [39]. Although not all bonds were formed again during the stabilization of the aggregates during the two wetting/drying cycles after the mortar treatment, we have samples with a better homogenization of the mucilage concentration, the porosity and a structural stability similar to undisturbed media.

The structural stability is more important in the presence of mucilage than in the absence of mucilage. Figure 4.9 illustrates this stronger structural stability in the presence of mucilage. In addition, we noted that the wetting time of the samples increased with the concentration of mucilage present in the medium. Moistening is immediate for concentrations of *c* = 0%. If *c* = 0.1%, humectation takes 2-3 seconds, whereas 3-5 minutes are need for *c* = 0.5% and an hour or even half a day for *c* = 1%. These observations confirm that mucilage concentration changes the hydraulic properties and structure of the medium.

<figure><img src="image_72.png" /><figcaption>Figure 4.9: Structural stability of samples at different concentrations from a medium modified after mixing substrate and mucilage.</figcaption></figure>

## 4.2 Characterization of hydraulic properties of the rhizosphere analogue medium

The aim of this characterization is to estimate the hydraulic conductivity and water retention of this analogue medium, with an emphasis on the wettability of the substrate in presence of amphiphilic (hydrophilic/hydrophobic).
**4.2.1 Materials and methods**

**4.2.1.1 Analogue rhizosphere sampling**

Several sample formats have been designed in order to perform the various experiments under optimum conditions for defining the hydraulic characteristics of the media. The design of each sample follows the steps described in section 4.1.3.

**Hydraulic conductivity** Samples from the disturbed medium design are used for Ksat measurements. In cylinders of 5 cm high and 5.1 cm in diameter, 100 g of dry sand with the appropriate concentration of dry mucilage is placed in the cylinder. The bottom of each cylinder is covered with muslin held in place with rubber bands, i.e the sample remains in place while the water drains away. We note that the dimensions of the cylinder require a large quantity of mucilage for this measurement. Then for these measurements, we selected two mucilage concentrations (0.1 and 0.5%) and a control medium without mucilage (0%). Three replicates were designed for the 0.1 and 0.5% concentrations and only two without mucilage. We have a total of 8 samples for Ksat measurements, whose height will be around 3.3 cm in the cylinders. Note that the cylinder must be filled at least to 3 cm of its height for the measurement to be carried out. Following the wetting/drying cycle, the dry samples are placed on a suction table to saturate them with water: they are now ready for measurement.

**Water retention** The disturbed medium design is used to create samples for measuring water retention. In small cylinders of 1 cm high and 2.6 cm in diameter, 7 g of sand with the desired mucilage concentration is placed in the cylinder. The bottom of the cylinder is closed with a cloth, allowing the sample to remain in place while the water drains away. We selected three mucilage concentrations (0.1%, 0.5% and 1%) in addition to control samples (0%). We performed two types of water retention measurements in these media. Initially, we made 10 replicates of each medium, and each replicate was used to measure water retention for the following suctions: 0.04, 0.11, 0.26, 0.51, 0.75, 0.99, 1, 2, 10 and 20 m. However, due to degradation issues of mucilage concentration over time, new samples were designed, i.e. 5 replicates of each medium for each suction to be measured: 0.25, 0.35, 0.5, 0.65, 1, 3 and 10 m. Thus, we have a first version with 40 samples designed while the second version requires 140 samples. The second version requires more samples and more mucilage, but limits mucilage degradation during measurements. Each sample will have a height of around 0.9 cm in the cylinders. Note that the cylinder must be filled to at least 5 mm of its height. Following the wetting/drying cycle, the dry samples are placed on a suction table to saturate them with water: they are then ready for measurement on the suction table and the pressure plate apparatus.

**Drop infiltration** Wettability measurements are performed by infiltrating a drop into two sample designs: disturbed and undisturbed medium. Small petri dishes of 1 cm high and 5.3 cm in diameter are used for disturbed media, while large petri dishes of 1 cm high and 8.7 cm in diameter are used for undisturbed media. Note that the bottom of each dish is closed by a rigid plastic base, allowing the sample to remain in place and retaining any water that infiltrates. We selected three mucilage concentrations (0.1%, 0.5% and 1%) in addition to control samples (0%) for each sample design. We have small petri dishes with 17 g of sand with mucilage concentration with the disturbed design and
large petri dishes with 30 g of sand with mucilage concentration with the undisturbed design. Note that each sample is approximately 0.52 cm in height. Following the wetting/drying cycle, the dry samples are ready for infiltration measurements.

### 4.2.1.2 Saturated soil hydraulic conductivity measurement

The hydraulic conductivity depends on the intrinsic properties of the porous medium, the degree of saturation and the properties of the flowing fluid. The presence of organic matter such as EPS from mucilage can alter the properties of the porous medium (pore distribution, grain wettability) as well as those of the fluid with a change in solutes (variation in viscosity, surface tension). As a result, the hydraulic conductivity is also altered. It is relevant to study the hydraulic conductivity of analogue

<figure><img src="image_73.png" /><figcaption>Figure 4.10: A constant head permeameter system.</figcaption></figure>

media for different mucilage concentrations.

The saturated hydraulic conductivity $K_{sat}$ is measure using a water system (Figure 4.10) This system is connected to the sample by pipes imposing a constant upward hydraulic head in the porous medium using a water height difference between two reservoirs. The water flowing out from the top of sample is collected on a precision balance placed under the evacuation tube to measure the accumulated mass of water and recorded each minute on a computer. Measurements at the same pressure continues until the flow rate is constant over several measurements. Several minutes can be necessary to reach equilibrium between 2 successive pressure steps. The water flow is measured at different pressures, and the $K_{sat}$ is estimated from linear regression based on Darcy law:

$$Q = K_{sat}S\frac{H}{L}, \qquad (4.1)$$

where $Q$ [g.min⁻¹] is the flow rate, $S$ (19.6 cm²) is the cross-section of the column, $K_{sat}$ [cm.min⁻¹] is saturated hydraulic conductivity, $H$ [cm] is the hydraulic head difference across the column, and $L$ (3.3 cm) is the height of the column.
4.2.1.3 Soil water retention curve

The water retention (i.e the soil moisture content for a given suction) of sand/mucilage media was measured for suctions varying in the range [0, 15] bars. Note that, the suction of 15 bars being the plant's permanent wilting point where it is considered that the water held in the soil at or below this suction is no more available for plant transpiration. According to the range of suction imposed to samples, two types of system were used: suction tables and pressure plate apparatus.

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

<figure><img src="image_75.png" /><figcaption>Figure 4.11: Water retention measurement system. (a) Suction table. (b) Pressure plate apparatus.</figcaption></figure>

**The suction table** In the low suction range from 0 to 0.1 bars, suction tables are used (Figure 4.11a). Samples are placed on a bed of saturated fine sand (silica) connected to a suction system. This suction is imposed by a differential water height between the water level in an adjustable water reservoir and the saturated surface of the sand bed. The sample are first saturated on the table, by rising the water level at the surface of sand or slightly above. Then negative heads are imposed by lowering the water reservoir level at a given height. The imposed head is maintained until hydraulic equilibrium is reached. In order to confirm equilibrium, a replicate of the samples is weighed every day until constant mass between 2 consecutive days. After equilibrium, the gravimetric water content of the sample (g of water per g of dry soil) is measured by oven drying 24 h at 105°C.

**The pressure plate apparatus** In the high-suction range from 0.1 to 15 bars, pressure plate device was used (Figure 4.11b). Samples are placed in a sealed enclosure on a micro-porous ceramic plate. The plate has been saturated and connected to the outside of the enclosure to evacuate any water excess. It is important that the ceramic plate and the soil sample are well hydraulically connected, i.e there is no air between samples and ceramic plate. When the chamber is closed, gas nitrogen is injected to impose a constant positive pressure. The pressure in the enclosure pushes part of the water contained in the initially saturated sample out of the enclosure. The water flow stops and soil water reaches equilibrium when the pressure (in absolute value) in the soil equals the gas pressure. In order to check the equilibrium of the sample, the evacuated water volume is estimated each day by recording the height of water in a small tube placed at chamber outlet. In addition, the tube can be weighed to
add a second control criterion. Equilibrium is considered to be reached when the mass of the tube and the level of water in the tube do not change for two consecutive days. When equilibrium is reached, samples are collected and their gravimetric water content measured.

#### 4.2.1.4 Wettability and drop infiltration measurement

The determination of the retention curve presented here is associated, in relation to wettability, with the hydrophilic properties of the porous medium. However, it is difficult to assess the hydrophobicity of a medium using this method. There are other standardized methods for measuring soil hydrophobicity. These include the Molarity of Ethanol Droplet (MED) test [357, 194] and the Water Drop Penetration Time (WDPT) test [357]. These methods classify the degree of hydrophobicity of the medium using criteria based on drop infiltration times. This measurement gives us a trend of the hydrophobicity in the medium.

Based on these tests, we aim to classify the wettability of our analogue media in order to characterize the effects of EPS concentration on the hydrophobicity level of the medium.

A water drop is deposited on our analogue media (disturbed and undisturbed) using the WDPT and MED tests to classify their degree of hydrophobicity (Figure 4.12).

<figure><img src="image_76.png" /><figcaption>Figure 4.12: Photo of droplet deposited on rhizosphere analogue media using an adjustable laboratory micropipette.</figcaption></figure>

**WDPT test** is a measure of hydrophobicity at the surface of a defined medium, based on the infiltration time of a drop into the medium deposited on its surface. A degree of hydrophobicity is associated with the drop's penetration time: instantaneous for a hydrophilic medium, very slow for hydrophobic media and staying at the surface until evaporation for media with very high hydrophobicity [338]. The definition of infiltration time associated with a hydrophobicity category is an arbitrary choice explaining the variety of classification tables present in the literature [32, 104, 103]. To estimate hydrophobicity in our analogue media (disturbed and undisturbed), we selected the classification defined by Bisdom *et al* (Table 4.2 for a 50 µL water drop). This classification is well described, and has equivalence with measurements of the contact angle of a drop on the medium. Contact angle is a common measure of surface wettability. This angle increases with hydrophobicity: an angle of 0° defines a hydrophilic medium, while an angle of 180° defines an extremely hydrophobic medium. However, it is difficult to
<table><caption>Table 4.2: Classification of the hydrophobicity of a medium: WDPT test.</caption><thead><tr><th></th><th>Bisdom et al., 1993 [32]</th><th>Mandal et al., 2009 [224]</th></tr><tr><th>Classification</th><th>Infiltration time [s])</th><th>Contact angle (°)</th></tr></thead><tbody><tr><td>Hydrophilic</td><td>&lt; 5</td><td>0</td></tr><tr><td>Slightly hydrophobic</td><td>5 à 60</td><td>67</td></tr><tr><td>Strongly hydrophobic</td><td>60 à 600</td><td>90</td></tr><tr><td>Severely hydrophobic</td><td>600 à 3600</td><td>98</td></tr><tr><td>Extremely hydrophobic</td><td>&gt; 3600</td><td>122</td></tr></tbody></table>

measure because it requires an equilibrium state which is difficult to obtain because of the dynamics of infiltration into the medium and its spreading over the surface.

MED test also known as the "Percentage Alcohol Test" [357, 194], is an indirect measure of soil surface tension. It also indicates the repulsion of a drop of water solution containing some percentage of ethanol from the soil surface. The presence of ethanol modifies the surface tension of the drop. The surface tension of the drop decreases as the percentage of ethanol in the solution increases. This test measures the volume percentage of ethanol in a drop of water required to reduce the drop's surface tension compared with the surface tension of the soil. The drop will rapidly penetrate the medium if its surface tension is lower than the surface tension of the medium. The spreading coefficient $S$ is then positive, i.e. we have a hydrophilic tendency promoting infiltration. Otherwise, the drop will penetrate the medium slowly. The spreading coefficient $S$ is then negative, i.e. we have a partial wetting tendency causing slower infiltration than with total wetting. The volume percentage of ethanol in a drop of water is associated with a drop penetration time. This time is compared with a reference infiltration time threshold $t_c$ to classify the degree of hydrophobicity at the surface of the medium. The reference time is an arbitrary criterion that varies in the literature: $t_c = 3$ s [83], $t_c = 5$ s [194] or $t_c = 10$ s [164].

The Crockford criterion has been chosen for our study in order to have classification equivalences between the MED test and the WDPT test. The penetration of a drop of 50 µL of water and ethanol solution is measured to determine the minimum volume percentage of ethanol that satisfies infiltration in less than 3 seconds. Each minimum volume percentage of ethanol is associated with a degree of hydrophobicity at the surface of the medium (Table 4.3 [103]). This degree of hydrophobicity is categorized into 7 classes: from 1 (very hydrophilic) to 7 (extremely hydrophobic).

<table><caption>Table 4.3: Classification of the hydrophobicity of a medium: MED test.</caption><thead><tr><th>Class</th><th>Description</th><th>Ethanol (%)</th></tr></thead><tbody><tr><td>1</td><td>Very hydrophilic</td><td>0</td></tr><tr><td>2</td><td>Hydrophilic</td><td>3</td></tr><tr><td>3</td><td>Slightly hydrophobic</td><td>5</td></tr><tr><td>4</td><td>Moderatly hydrophobic</td><td>8.5</td></tr><tr><td>5</td><td>Strongly hydrophobic</td><td>13</td></tr><tr><td>6</td><td>Very strongly hydrophobic</td><td>24</td></tr><tr><td>7</td><td>Extremely hydrophobic</td><td>36</td></tr></tbody></table>
4.2.2 Results

4.2.2.1 Saturated hydraulic conductivity

Saturated hydraulic conductivity is measured from average of two samples for the control medium (c = 0%) and of three samples for the other two concentrations (c = 0.1% and c = 0.5%).

Without mucilage (c = 0%), the measured Ksat is $1.35 \times 10^{-4}$ m.s$^{-1}$, and quite close to mea-
surements obtained from independent experiments in the literature (Table 4.1). For higher mucilage
concentrations, we found that Ksat was of the same order of magnitude as the control medium. More-
over, there is a slight increasing tendency of Ksat with increasing mucilage concentration: Ksat is
$1.94 \times 10^{-4}$ m.s$^{-1}$ for c = 0.1% and Ksat is $2.26 \times 10^{-4}$ m.s$^{-1}$ for c = 0.5%. Our results suggest
that the presence of molecules adsorbed on the particles and/or diffusing in the liquid did not affect
hydraulic conductivity.

Our results in the presence of mucilage are surprising compared with the trend obtained in the
work of Kroener et al [199]. They showed that mucilage influences the hydraulic conductivity of
the medium as a function of particle size. In their experiment, they produced 5 ranges of analogue
rhizosphere medium from 5 different ranges of particle sizes, from coarse sand to silt (630 – 1000,
360 – 630, 200 – 500, 60 – 200, and < 20 µm) were prepared with varying concentrations of chia
seed mucilage (0%, 0.1%, 0.4%, and 0.8%). Note that their medium design seems to be similar to
our undisturbed design, except for the medium with particle sizes < 20 µm. This latter medium was
sieved to break aggregates to enable homogeneous mixing with mucilage. This process is similar to
disturbed design of our samples. It should also be noted that the substrates chosen in their experiment
were not Fontainebleau sand. They made 3 replicas of each medium: two media using quartz sand
from sand pit in southern Florida, USA (coarse sand with a particle grain size of 630 – 1000 µm and
360 – 630 µm), two other media using quartz sand from sand pit located near Duingen, Germany
(fine sand with a particle grain size of 200 – 500 µm and 60 – 200 µm) and another using silt from
Reinhausen, Lower Saxony, Germany (particle grain size < 20 µm). The Fontainebleau sand (fine sand
with a particle grain size of 200 – 250 µm) and the sand near Duingen have similar particle grain sizes.
We can therefore expect similar results for this range of media.

However, they obtained a decreasing trend in hydraulic conductivity at saturation for all sand
media with increasing mucilage concentration. Note that this trend is similar with the results in Landl
et al [206]. For every experiment, their hydraulic conductivities for sandy media decreased by several
orders of magnitude:

* Ksat was $3 \times 10^{-4}$ m.s$^{-1}$ for c = 0%, whereas it was $5 \times 10^{-7}$ m.s$^{-1}$ for c = 0.8% in the medium with a particle size 200 – 500 µm in the results of Kroener,

* Ksat was $10^{-4}$ m.s$^{-1}$ for c = 0%, whereas it was around $10^{-6}$ m.s$^{-1}$ in presence of EPS in the results of Landl.

Nevertheless, we note that the silt medium, in the Kroener experiment, which seems to have a design
close to that of our media, has a similar trend to our results. The hydraulic conductivity at saturation
remains in the same order of magnitude as the medium without mucilage (Figure 4.13a). Comparing
the ratio between the saturation conductivity for a EPS concentration c and the saturation conductivity
without EPS, we observe that the relative hydraulic conductivity of the Fontainebleau sand (orange
dotted line) is similar to that of the silty medium (black solid line) in the experiment of Kroener (Figure 4.13b).

We observe that, depending on the design of the analogue rhizosphere medium, results on hydraulic conductivity can vary, as shown by our measurements on Fontainebleau sand and sand near Duingen, which have the same particle size but not necessary the same bulk density. Due to this variability, we will assume in our characterization that mucilage does not affect the saturation hydraulic conductivity of our medium. Further measurements are required to confirm and quantify the effect of mucilage on hydraulic conductivity. Moreover, it should be interesting to compare the bulk density of each medium that can affect the $K_{sat}$ value.

<figure><img src="image_77.png" /><figcaption>Figure 4.13: Saturated hydraulic conductivity (a) as a function of mucilage concentration for various ranges of soil particle sizes, and (b) the relative reduction of saturated hydraulic conductivity as a function of mucilage concentration for various ranges of soil particle sizes. Solid lines correspond to Kroener *et al* data [199] and dotted lines correspond to our measurements.</figcaption></figure>

### 4.2.2.2 Retention curve following single samples

In this experiment 1, we selected suctions between 0 and 20 m. At the end of these measurements, the dry density of our analogue media was similar for the different concentrations ($\gamma_d$ = 1.31). The results show that porosity is similar in all designed samples. The porosity $\phi$ = 0.33 is more or less the same as that obtained in the literature ($\phi$ = 0.36). Figure 4.14 represents the average retention curves obtained for our analogue media with $c$ = 0, 0.1, 0.5 and 1%. This figure shows water content results expressed as gravimetric water content (Figure 4.14a) or water saturation (Figure 4.14b).

We observe trends in retention depending on the imposed suctions and EPS concentrations. For a given EPS concentration, moisture content in the medium follows the classical trend of a decrease with increasing suction. However, for a given suction, the water content increases with concentration. The attractive forces of the EPS in the medium enable more water to be absorbed in the sandy soil and a higher retention. This trend is visible all over the range of investigated suctions, except for the suction between 0.26 and 0.75 m, where this trend is reversed. As EPS concentration increases, the medium retains less water for these suctions. This could be interpreted as a greater repulsion of water by the EPS, with a greater proportion of repellent EPS parts compared to attractive parts, in this zone of sharp decrease in water content. So EPS could induce 2 kind of trends for retention: higher retention for low and high suctions (hydrophilic behaviour) and a lower retention for intermediate
<figure><img src="image_78.png" /><figcaption>Figure 4.14: Average retention curve of our analogue media for concentrations of $c = 0, 0.1, 0.5$ and $1$, using the non-destructive method: (a) gravimetric water content and (b) saturation.</figcaption></figure>

suctions (hydrophobic behaviour), but with highly varying water content.

This trend is consistent with the concept of a non-monotonic retention curve, difficult to experimentally determine due to the processes involved in measuring the curve and multiple equilibrium states. The retention curve is not unique. It varies according to the measurement carried out during humidification or drainage (as in this case). The hysteresis nature of the retention curve confirms that it is not monotonic, particularly in the presence of amphiphilic matter.

<figure><img src="image_79.png" /><figcaption>Figure 4.15: Photo of a Collembola on our sample support after passing through the pressure plate apparatus.</figcaption></figure>

However, this result should be linked to the duration of the experiment (5 months) to obtain the hydraulic equilibrium for each suctions ranging from 0 to 20 m. During this period, the mucilage was probably modified due to environmental conditions. Indeed, it is impossible to maintain a sterile experimental environment throughout this period in our lab conditions. This observation is supported by the presence of Collembola in our samples (Figure 4.15). The presence of these organisms question about time evolution of EPS concentrations in the medium.

The origin of the Collembola development could be attributed to the soil used to enhance a hydraulic connection of the samples to the porous plate in the pressure plate apparatus. This soil must certainly
contain organisms. In this case, it is important to ensure that the change in retention curve trend is
due to hydrophobicity and not the presence of the organisms degrading the mucilage. We propose to
reproduce the previous experiment, adding more control over the sterility of the system and carrying
out measurements on a single sample for each pressure studied. In addition, we are adding a time
restriction on water equilibrium in the sample to limit mucilage degradation.

4.2.2.3 Retention curve with the "one pressure - one sample" process

In order to limit degradation of mucilage in the medium, instead of following a single sample through
multiple pressure steps necessitating long times, we used one sample for one pressure step. The time for
measurement is considerably reduced, reducing degradation impacts. Globally, hydraulic equilibrium
of samples was reached in less than 12 days.

In parallel, measurements of microbial respiration and the amount of carbon present in the media
were determined over time in the mucilage sand media. Results show that microbial activity was
present in our media, causing the mucilage concentration to vary with time (Appendix G). A peak
in respiration activity was detected at 12 days. This experiment shows that the biological activity
complexifies the medium characterization, enhancing the need for a balance of time between time for
sample equilibrium and degradation.

<figure><img src="image_80.png" /><figcaption>Figure 4.16: Average retention curve of our analogue media for concentrations of c = 0, 0.1, 0.5 and 1, using the destructive method: (a) gravimetric water content and (b) saturation.</figcaption></figure>

In this second water retention experiment, we selected pressures mainly for low suctions, to confirm
the inversion in mucilage effect shown above.

Figure 4.16 represents the average retention curves obtained for our analogue media for mucilage concentrations of c = 0, 0.1, 0.5 and 1%.

The trends obtained in the previous determined retention curve are again present, with inversion
of mucilage effect for suctions between 0.35 and 0.65 m. This result confirms the idea that repulsive
actions act in the presence of amphiphilic matter at an intermediate level of saturation.

4.2.2.4 Wettability measurement and drop infiltration test

WDPT test

We observe that drop infiltration time increases with EPS concentration in both types of media. Media become more hydrophobic in the presence of EPS. Comparing the level of hydrophobicity results for the same concentration, we observe high differences between the two types of samples. Infiltration takes more time for the undisturbed medium. Figure 4.17 displays the average infiltration time of a water drop in our dry analogue media (disturbed and undisturbed samples for concentrations of c = 0, c = 0.5 and c = 1%). The undisturbed medium at c = 0.5% is severely hydrophobic (t ≈ 3500 s) and at c = 1% is on the boundary between severely and extremely hydrophobic (t ≈ 3600 s). In contrast, the disturbed medium at c = 0.5% is hydrophilic (t ≈ 2 s) and at c = 1% is strongly hydrophobic (t ≈ 380 s). Undisturbed media appear more hydrophobic than disturbed media with the same EPS concentration. This difference seems due to a difference in distribution of EPS in the undisturbed media. A crust is present on the surface of these undisturbed samples, in contrast to disturbed. The crust results from the accumulation of EPS on the surface of the media, making it more hydrophobic.

<figure><img src="image_81.png" /><figcaption>Figure 4.17: Infiltration time of a drop in disturbed sample (3 replications) compared with the infil-
tration time in a undisturbed sample (2 replications).</figcaption></figure>

If the WDPT test is a suitable and reproducible method for our laboratory measurements it can be
relatively long especially for highly hydrophobic media (e.g. for c = 0.5% and c = 1% in undisturbed
media). In this case, evaporation can affect our infiltration time measurement. Despite this constraint,
we have identified a trend of increasing hydrophobicity at the surface of our media with increasing
EPS concentration.
MED test

Drop infiltrations with different percentages of ethanol by volume were carried out in our undisturbed samples exhibiting the high hydrophobicity observed with the WDPT test (c = 0.5% and c = 1%). The c = 0.5% and c = 1% media were found to be consistent with the classification defined by the WDPT test: the c = 0.5% medium is classified as a very strongly hydrophobic medium (24% ethanol) and the c = 1% medium is classified as an extremely hydrophobic medium (36% ethanol). The MED test confirms the trend obtained with the WDPT test: hydrophobicity increases with EPS concentration.

This MED test is more suitable for highly hydrophobic media than WDPT test. It gives a rapid result, limiting the effects of evaporation, especially in field measurements where it is difficult to control external phenomena impacting the infiltration of a drop [83, 352, 103]. However, in our situation, the WDPT test remains the most practical and reproducible because our measurements are carried out in a laboratory (controlled environment) [83, 103] and there is no addition of any other substance than water to the medium. This last point is particularly important for comparing the infiltration time between the successive infiltration of one and two drops, as addition of second drop occurs on an already moistened medium.

Hydrophobicity and water saturation in the presence of EPS from successive drop infil-
tration test

Hydrophobicity increases with the presence of EPS at the surface of dry medium. We studied the infiltration behavior of a second drop placed at the same location as the first one in the disturbed medium. Figure 4.18 shows the infiltration time of the first and then second drop of the same volume in our disturbed media with concentration c = 0%, c = 0.1%, c = 0.5% and c = 1%. We observe that

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

Figure 4.18: Infiltration time of a first drop in a disturbed sample (3 replications) compared with the infiltration time of a second drop in the same sample. * Student's *t* test (paired samples) shows a significant difference (p<0.0001) between the infiltration time of a second drop compared with the first drop, when the concentration of EPS is 1% in the sample tested.
the second drop needs more time to infiltrate than the first. Apart the lowering of suction gradients in the second drop experiment this also means that the surface of the medium becomes more hydrophobic after the first drop. This is particularly illustrated with the 1% concentration where the first drop infiltrates the medium in around 380 s compared to 1130 s for the second. This increase in the time infiltration of the second drop is similar for each concentration compared with the infiltration time of the first drop on dry soil. Moreover, this infiltration time increases gradually with EPS concentration for both the first and second drops.

Hydrophobic interactions between the EPS and the droplet are always present that increase time infiltration: lateral and vertical diffusion in the medium is subject to these hydrophobic interactions. At the surface, the second drop spreads more than the first one due to the effect of attractive interactions. But the majority of EPS in soil are still hydrophobic at low saturation [105], which increases infiltration time. Thus the degree of hydrophobicity changes with the water saturation at the surface of the medium in the presence of EPS. Wettability of the medium varies with saturation at the surface of the medium [105, 219]. Hydrophobicity increases after drying a very wet sample [230, 134]. The wettability of our media in the presence of EPS changes depending on their concentration and initial saturation in the medium.

## 4.3 Parametrization of soil hydraulic properties accounting for variable wettability

Modeling of water transfer in porous media requires relationships which describe the hydraulic properties of the medium. These relationships are defined by the kind of (semi-empirical) function chosen and its associated parameters (chapters 2-3). The aim of this section is to define the parameters of the hydraulic functions associated with our disturbed media. These parameters are fitted based on empirical data obtained during the characterization of hydraulic properties (section 4.2).

### 4.3.1 Wetting function: disjoining pressure

Wettability in the analogue medium is modeled by the disjoining pressure $\Pi_w$. This pressure is the linear combination of the attractive $h_c$, repulsive $h_r$ and compressive $h_{comp}$ pressures (Chapter 2),

$$ \Pi_w(s) = h_c(s) + h_r(s) + h_{comp}(s). \qquad (4.2) $$

The retention curve during drying is associated with attractive pressure (section 4.2.2.2). This curve is generally presented by the Brooks and Corey function [42]

$$ h_c(s) = H_{cap}s^{-\frac{1}{\lambda}}, \qquad (4.3) $$

or the van Genuchten-Mualem function [240, 343]

$$ h_c(s) = \frac{1}{\alpha} \left( s^{-\frac{1}{m}} - 1 \right)^{\frac{1}{n}}, \qquad (4.4) $$

with $m = 1 - \frac{1}{n}$.
These two functions have two parameters that can be fitted to data using the least squares method. The fitted parameters are given in Table 4.4. Figure 4.19 compares the experimental retention results

<table><caption>Table 4.4: Parameters Fitted to the characterization of soil water properties in presence of amphiphilic matter. The symbol "X" corresponds that the data is not defined by characterization.</caption><thead><tr><th>Experiment n°</th><th>c</th><th>λ</th><th>φ</th><th>H<sub>cap</sub> [m]</th><th>K<sub>sat</sub> [m.s<sup>-1</sup>]</th><th>σ</th><th>s<sub>r</sub></th><th>κ</th><th>m</th><th>n</th><th>α [m<sup>-1</sup>]</th></tr></thead><tbody><tr><td>Exp 1</td><td>0</td><td>0.87</td><td>0.33</td><td>0.118</td><td>1.35 × 10<sup>-4</sup></td><td>X</td><td>X</td><td>X</td><td>0.833</td><td>5.98</td><td>2.24</td></tr><tr><td>Exp 1</td><td>0.1</td><td>0.83</td><td>0.33</td><td>0.102</td><td>1.94 × 10<sup>-4</sup></td><td>X</td><td>X</td><td>X</td><td>0.8</td><td>5</td><td>2.4</td></tr><tr><td>Exp 1</td><td>0.5</td><td>1.01</td><td>0.34</td><td>0.14</td><td>2.26 × 10<sup>-4</sup></td><td>X</td><td>X</td><td>X</td><td>0.8</td><td>5</td><td>2.49</td></tr><tr><td>Exp 1</td><td>1</td><td>0.74</td><td>0.36</td><td>0.104</td><td>X</td><td>X</td><td>X</td><td>X</td><td>0.8</td><td>5</td><td>2.81</td></tr><tr><td>Exp 2</td><td>0</td><td>1.23</td><td>0.52</td><td>0.24</td><td>X</td><td>X</td><td>X</td><td>X</td><td>0.685</td><td>3.16</td><td>2.11</td></tr><tr><td>Exp 2</td><td>0.1</td><td>0.86</td><td>0.55</td><td>0.16</td><td>X</td><td>X</td><td>X</td><td>X</td><td>0.685</td><td>3.16</td><td>2.11</td></tr><tr><td>Exp 2</td><td>0.5</td><td>0.51</td><td>0.69</td><td>0.046</td><td>X</td><td>X</td><td>X</td><td>X</td><td>0.686</td><td>3.19</td><td>4.17</td></tr><tr><td>Exp 2</td><td>1</td><td>0.75</td><td>0.43</td><td>0.182</td><td>X</td><td>X</td><td>X</td><td>X</td><td>0.698</td><td>3.31</td><td>2.65</td></tr></tbody></table>

for each concentration obtained for the measurement following single samples (Exp 1 in section 4.2.2.2) with the parameter fits for the Brooks and Corey (BC) and van Genuchten (VG) function. Figure 4.20

<figure><img src="image_83.png" /><figcaption>Figure 4.19: Measured and fitted water retention curve for sandy soil with and without EPS concentration using Brooks and Corey (BC) and Van Genuchten (VG) functions. Measurements as well as fitted lines according to the parameters of Table 4.4 are shown in Exp 1.</figcaption></figure>

compares the experimental retention results for each concentration obtained for the measurement using "one pressure - one sample" process (Exp 2 in section 4.2.2.3) with the same fitting function. The fit with the Brooks and Corey model is found to be better adapted to the experimental data for each concentration and we will use it with the parameters determined in Table 4.4 for the description of the attractive pressure.
<figure><img src="image_84.png" /><figcaption>Figure 4.20: Measured and fitted water retention curve for sandy soil with and without EPS concentration using Brooks and Corey (BC) and Van Genuchten (VG) functions. Measurements as well as fitted lines according to the parameters of Table 4.4 are shown in Exp 2.</figcaption></figure>

The repulsive pressure is based on a function of embedding 4 parameters, describing its influence in the intermediate saturation range. Based on the changes observed in the retention curve, we consider that this repulsive pressure acts for suctions ranging from 0.35 to 0.65 m (section 4.2.2.3) corresponding to saturation between $s = 0.2$ and $s = 0.6$. We represent the repulsive pressure with a Gaussian function with a standard deviation of $\sigma = 0.2$ and residual saturation given by $s_r = 0.05$ then a maximum repulsion occurs around $s = 0.2$.

$$h_r(s) = -\frac{\chi}{H_{cap}\sigma\sqrt{2\pi}} \exp\left(-\frac{(s - s_r)^2}{2\sigma^2}\right), \qquad (4.5)$$

where $\chi$ is the amplitude of repulsion.

In the case of sand, research indicates hydrophobicity effects for saturation around $s = 0.2$ [219]. This result is consistent with our function. However, the parameters determined here are based on tendencies and need to be defined more accurately with appropriate experiments, in particular for the amplitude of repulsion. This parameter should vary according to the concentration of EPS in the medium. A capillary rise experiment in a stratified medium (one stratum without EPS and another with different EPS concentrations) [201] could enable us to define the $\chi$ parameter for our media.

Finally, we define an air compression function $h_{comp}$ taking into account air trapping when the medium is close to saturation. This function is defined by Cueto and Juanes [86] by the following function:

$$h_{comp}(s) = \operatorname{sgn}(1 - \lambda) H_{cap} \exp(-\kappa(1 - s)) s^{-\frac{1}{\lambda}} \left(1 + \frac{\lambda}{\lambda - 1} \kappa s\right). \qquad (4.6)$$
It depends on two parameters: a coefficient $\kappa$ and the shape parameter $\lambda$ derived from the retention curve. The $\kappa$ coefficient is associated to the maximum saturation considered in our study, i.e. we define $\kappa$ such that we have a maximum water saturation in our sample $s_{comp}$ close to 1 ($s_{comp} \approx 1$). This model allows us to take account of air trapping in the medium.

We have determined most of the wettability parameters from our characterizations of disturbed analogue rhizosphere media. The remaining parameters $\chi$ and $\kappa$ according to the EPS concentration need to be characterized and parameterized. Literature experiments will be used to complete these parameters required for numerical simulations (Chapter 5).

### 4.3.2 Hydraulic conductivity

Hydraulic conductivity in the disturbed medium is modeled using Eq. (1.27) with a relative conductivity function $k_r$ and the hydraulic conductivity at saturation $K_{sat}$. Hydraulic conductivity at saturation is obtained from measurements in section 4.2.2.1. The relative conductivity is represented by the van Genuchten-Mualem function [343, 240],

$$k_r(s) = \sqrt{s} \left[ 1 - \left( 1 - s^{\frac{1}{m}} \right)^m \right]^2, \quad (4.7)$$

with $m = 1 - \frac{1}{n}$.

This function can be calibrated more accurately by measuring relative conductivity. However, we believe that this measurement is not the most important, because the effects of relative conductivity are less important than the pressures associated with the wettability of the medium. Parameters of the relative conductivity function depend on the shape parameter $n$ ($m = 1 - \frac{1}{n}$) of the drying water retention curve Eq. (4.4) that was determined using Van Genuchten parameter defined in Table 4.4. Moreover, the relative conductivities for the different concentrations are more or less the same (Figure 4.21).

<figure><img src="image_85.png" /><figcaption>Figure 4.21: Relative conductivity of (a) undisturbed and (b) disturbed media for the concentration $c = 0\%$, $c = 0.1\%$, $c = 0.5\%$ and $c = 1\%$.</figcaption></figure>

Therefore, the hydraulic conductivity function is described according to the van Genuchten-Mualem function with the measured conductivity at saturation and the parameters of the retention curve.
4.4 Conclusion

We presented the creation of rhizosphere analogue media from Fontainebleau sand and Chia seed mu-
cilage for EPS concentration ranges referenced in a rhizosphere soil (c = 0%, c = 0.1%, c = 0.5 %
and c = 1%). On these media, we carried out measurements to characterize its hydraulic proper-
ties: measuring hydraulic conductivity at saturation, defining the drying water retention curve and
defining drop infiltration trends. During these measurements, we noted that temporal constraints due
to the equilibrium of the water content in the medium, as well as the time degradation of EPS by
microbial activity, had to be taken into account (Appendix G). Moreover, we highlight the fact that
mucilage quantities are generally rather low, which limits the number of replicates in the statistical
study. Despite this, we noted that hydrophobicity in the media increased with EPS concentration and
influenced the drying water retention curve: variable trend observed as a function of saturation and
EPS concentration in the medium.

Nevertheless, we did not observe any significant variation in hydraulic conductivity at saturation, unlike Kroener *et al.* [199] experiments. We suspect that the design of the medium has a strong influence on its characterization, as we observed between disturbed and undisturbed soils. Finally, we have defined hydraulic parameters for the numerical simulation of analogue rhizosphere disturbed soil and different EPS concentrations using these experiments (Table 4.4).

Some parameters remain to be defined more precisely, in particular those related to the hydropho-
bicity of the medium. Additional experiments are required to parameterize them as capillary rise in a
stratified medium to parameterize the amplitude of repulsion, or infiltration in a stratified medium to
detect possible typical sizes of fingerings.

The water retention curves parameters obtained in this chapter are similar to those defined in Kroener *et al.* [201, 199]. However, they defined an EPS concentration thresholds illustrating notions of greater or lesser hydrophobicity (see section 5.2). In the next chapter, we will use these results to calibrate our disturbed media for simulating transfer in rhizosphere analogue media (infiltration of a drop and water uptake by a single root).
# Chapter 5

## Water transfer model in the rhizosphere in the presence of EPS

The aim of this chapter is to model water transfer in a rhizosphere analogue soil containing EPS. We propose to study two situations using two different models. These models are based on the same free energy representation of the water dynamics in the soil. The first model represents the dynamic of a drop infiltrating in a soil with different EPS concentrations. It takes into account several processes: the spreading of the drop on the surface of the soil, the infiltration and diffusion of the drop into the soil. The second model simulates water transfer in soil in the presence of a single root. The simulations with this model highlight the delays rewetting of the rhizosphere [366] compared to the bulk soil.

The first model coupling drop spreading and infiltration on the soil surface and drop in the soil is presented in section 5.1. The repellent energy density used in this model is calibrated depending on different EPS concentrations in the section 5.2. The coexisting equilibrium states of the saturation and the height of the drop are presented in section 5.3. This results obtained by this first model are discussed in section 5.4. The second model simulating water transfer in the soil in the presence of a single root is presented in section 5.5. Simulation of water transfer from bulk soil to the root obtained by this model are presented in section 5.6.

### Contents

<table><tr><td><strong>5.1</strong></td><td><strong>Model of drop infiltration in rhizosphere analogue soil</strong></td><td><strong>130</strong></td></tr><tr><td>5.1.1</td><td>Free energy of the drop/soil system</td><td>130</td></tr><tr><td>5.1.2</td><td>Drop on an impervious substrate</td><td>131</td></tr><tr><td>5.1.3</td><td>Unsaturated porous medium</td><td>133</td></tr><tr><td>5.1.4</td><td>Coupled drop/soil system</td><td>135</td></tr><tr><td>5.1.4.1</td><td>Interface Soil/Drop</td><td>135</td></tr><tr><td>5.1.4.2</td><td>Exchange of water</td><td>137</td></tr><tr><td>5.1.4.3</td><td>Governing equations of the coupled drop/soil system</td><td>137</td></tr><tr><td><strong>5.2</strong></td><td><strong>Calibration of soil model</strong></td><td><strong>138</strong></td></tr><tr><td><strong>5.3</strong></td><td><strong>Equilibria of flat film and uniform saturation</strong></td><td><strong>140</strong></td></tr><tr><td>5.3.1</td><td>Equilibria for fixed H</td><td>141</td></tr><tr><td>5.3.2</td><td>Equilibrium with respect to H</td><td>142</td></tr></table>
<table><tbody><tr><td><strong>5.4</strong></td><td><strong>Drop infiltration</strong></td><td><strong>143</strong></td></tr><tr><td>5.4.1</td><td>Wet porous medium</td><td>143</td></tr><tr><td>5.4.2</td><td>Dry porous medium</td><td>144</td></tr><tr><td>5.4.3</td><td>Pre-wetted porous medium</td><td>146</td></tr><tr><td>5.4.4</td><td>Conclusion</td><td>148</td></tr><tr><td><strong>5.5</strong></td><td><strong>Water transfer around a single root</strong></td><td><strong>148</strong></td></tr><tr><td>5.5.1</td><td>Model of soil around a single root</td><td>148</td></tr><tr><td>5.5.2</td><td>Model of water transfer around a single root</td><td>150</td></tr><tr><td><strong>5.6</strong></td><td><strong>Humectation of the soil</strong></td><td><strong>151</strong></td></tr><tr><td><strong>5.7</strong></td><td><strong>Conclusion</strong></td><td><strong>154</strong></td></tr></tbody></table>

## 5.1 Model of drop infiltration in rhizosphere analogue soil

We consider a planar 2D drop on a porous medium with a concentration $c$ of amphiphilic matter. The drop is defined by a water height $\zeta(x, t)$ and the porous medium contains a mean saturation profile $s(x, t)$ over a depth $D$. We define an equivalent height of water $\xi(x, t)$ corresponding to the mean saturation to simply the study. A sketch of the considered geometry is shown in Figure 5.1.

<figure><img src="image_86.png" /><figcaption>Figure 5.1: Sketch of the considered geometry for a liquid drop on a soil with amphiphilic matter.</figcaption></figure>

This situation corresponds to the wettability infiltration test used in the chapter 4. We could also related this situation to the rain drop infiltration into the soil.

### 5.1.1 Free energy of the drop/soil system

In this study, we neglect the evaporation in the different phases. We assume an Helmholtz free energy $\mathcal{F}$ in the closed system. The free energy is the sum of three terms: free energy in the film $\mathcal{F}^{\zeta}$, in the soil $\mathcal{F}^{\xi}$ and at the interface film/soil $\mathcal{F}^{\text{int}}$:

$$ \mathcal{F} = \mathcal{F}^{\zeta} + \mathcal{F}^{\xi} + \mathcal{F}^{\text{int}}. \qquad (5.1) $$

Each energy term $\mathcal{F}^{\square}$ has an associated $f^{\square}$ energy density and then

$$ \mathcal{F} = \int_{\Omega} f \, d\Omega = \int_{\Omega} [f^{\zeta} + f^{\xi} + f^{\text{int}}] \, d\Omega \qquad (5.2) $$

We called *wetting energy* of liquid film, resp. water in porous media, the energy resulting from the liquid/substrate interaction, resp. liquid/solid matrix interaction. The wetting energy is distinguished
from the free surface energy related to the surface tension. In fact, the latter is non-local since it depends on the derivative of ζ or ξ. Therefore, when we neglect gravity energy, the energy density of the system is given by

$$f = f_w + f_{\nabla}. \qquad (5.3)$$

Note that $f_w$ and $f_{\nabla}$ are also components of the three energy, then we note $f_w^{\square}$ and $f_{\nabla}^{\square}$ where $\square$ denotes the corresponding energy.

In the next subsections, we describe the model in each phase: drop on a substrate (Section 5.1.2) and water in porous medium (Section 5.1.3) separately. Then, we present the coupled system (Section 5.1.4) with respect to the free energy of the system.

## 5.1.2 Drop on an impervious substrate

In this section, we present the model of a thin film liquid on a flat impervious solid substrate in the case of partial wetting. A similar model is derived for the soil in the next section. The dynamics of the thin film is based on lubrication-type model, i.e. the flow is assumed be quasi-static approximation with long wavelength hypothesis.

Thin film model is based on mass conservation of the water height $ζ(\mathbf{x}, t)$:

$$\partial_t \zeta = -\nabla \cdot \mathcal{J}^{\zeta} \qquad (5.4)$$

where the flux $\mathcal{J}^{\zeta}$ is proportional to the pressure gradient $\nabla P$: $\mathcal{J} = Q(-\nabla P)$. Mobility $Q$, for no slip boundary condition at liquid/solid interface, depends on $ζ$ as follows:

$$Q(\zeta) = \frac{\zeta^3}{3\mu}. \qquad (5.5)$$

The pressure is the sum of three terms: the hydro-static pressure, the Laplace pressure, and the Derjaguin or disjoining pressure representing the effective molecular interaction between the substrate and the free surface [98, 93]. In the context of longwave approximation, the Laplace pressure is linearized and then

$$-P = \gamma\Delta\zeta + \Pi(\zeta). \qquad (5.6)$$

The parameters $ρ, g, \mu$ and $γ$ are the liquid density, the gravity, the dynamic viscosity and the liquid-gas interface tension, respectively. For submillimetric droplets, hydro-static pressure is negligible compared to the Laplace pressure or disjoining pressure. Here, we have decided to neglect hydro-static pressure, as its effect does not significantly affect drop behavior. From Eqs. (5.4) and (5.6), the governing equation of the film height can be written [257, 236]

$$\partial_t \zeta + \nabla \cdot [Q(\zeta)\nabla(\gamma\Delta\zeta + \Pi(\zeta))] = 0. \qquad (5.7)$$

This equation is known as the lubrication equation [22, 21] and it is proven to reproduce thin film dynamics for partial wetting. Partial wetting results from a competition between attractive and repulsive interaction acting at different scales between the substrate and the free surface [328, 18]. Usually, the
disjoining pressure is associated to the coexistence of an ultra-thin film $\zeta_p$ with a thick film defined by $\Pi(\zeta_p) = 0$ and the stability condition $\Pi'(\zeta_p) < 0$. This ultra-thin refers initially to a precursor film of a few nanometers which is very quickly developed before a drop completely spread. In this modeling the precursor height $\zeta_p$ represents a cut-off scale and not necessary a pseudo partial wetting. For instance, the static contact angle is related to $\Pi$ by the augmented Young-Laplace equation:

$$ \int_{\zeta_p}^{+\infty} \Pi(\zeta) d\zeta = \gamma (\cos \theta - 1). \qquad (5.8) $$

Among the expressions of the disjoining pressure with precursor [98, 326, 327, 22, 323, 155], we employ the following one [327, 22]:

$$ \Pi(\zeta) = \Pi_0 \left[ b \left( \frac{\zeta_p}{\zeta} \right)^3 - \exp \left( -\frac{\zeta}{\zeta_p} \right) \right] \qquad (5.9) $$

where a long range attractive interaction as $\zeta^{-3}$ and a short range exponential repulsive interaction. Note at very short range, the interaction is again attractive allowing the existence of stable precursor film. Three parameters characterize this disjoining pressure: $\Pi_0$ the pressure amplitude, $b$ dimensionless amplitude of the repellent interaction and $\zeta_p$ is the precursor film height. Depending on these

<figure><img src="image_87.png" /><figcaption>Figure 5.2: Wetting pressure in the thin film. Parameters are $b = 0.3679$, $\zeta_p = 0.02$ mm and $\Pi_0 = 9514$ Pa.</figcaption></figure>

parameters, the disjoining pressure $\Pi(\zeta)$ may have two zeros as displayed in Figure 5.2 $\zeta_1 \approx 0.02$ mm and $\zeta_2 \approx 0.134$ mm. The smallest corresponds to the precursor height $\zeta_1 = \zeta_p$. The interval $[\zeta_1, \zeta_2]$ defines the hydrophobic range, i.e. a thin film of height $\zeta \in [\zeta_1, \zeta_2]$ the wetting energy is minimized by decreasing the film height. Indeed, the disjoining pressure $\Pi$ is related to the wetting energy $f_w^\zeta$ in [92] by:

$$ \partial_\zeta f_w^\zeta(\zeta) = -\Pi(\zeta), \qquad (5.10) $$

where the wetting energy $f_w^\zeta$ can be express as the composition of attractive $f_{wa}^\zeta$ and repellent $f_{wr}^\zeta$
energies:

$$f_w^\zeta(\zeta) = f_{wa}^\zeta(\zeta) + f_{wr}^\zeta(\zeta), \qquad (5.11a)$$

$$f_{wa}^{\zeta}(\zeta) = \frac{\Pi_0 b \zeta_p^3}{2 \zeta^2}, \qquad (5.11b)$$

$$f_{wr}^{\zeta}(\zeta) = -\Pi_0 \zeta_p \exp\left(-\frac{\zeta}{\zeta_p}\right). \qquad (5.11c)$$

Eq.(5.7) can be written as a gradient dynamics of a conserved order parameter field [236, 271, 321]:

$$\partial_t \zeta = \nabla \cdot \left[ Q(\zeta) \nabla \frac{\delta \mathcal{F}^\zeta}{\delta \zeta} \right] \qquad (5.12a)$$

$$\mathcal{F}^\zeta[\zeta] = \int \left[ \frac{\gamma}{2} |\nabla \zeta|^2 + f_w^\zeta(\zeta) \right] d^2x, \qquad (5.12b)$$

$$\frac{\delta \mathcal{F}^\zeta}{\delta \zeta} = \partial_\zeta f_w^\zeta(\zeta) - \gamma \Delta \zeta. \qquad (5.12c)$$

where $\frac{\delta}{\delta\zeta}$ denotes the variational derivative and $\mathcal{F}^\zeta$ is the thin film free energy in long-wave approximation [375, 242]. The associated energy density contains the energy of the free surface (first term in Eq. (5.12b)) and the wetting energy (second term in Eq. (5.12b)). The variational derivative $\frac{\delta \mathcal{F}^\zeta}{\delta \zeta}$ is the pressure $P$ in the liquid: the Laplace pressure and the disjoining pressure defined by Eq (5.10).

This representation based on free energy represents dynamics associated with thin films on a substrate. In the following, we will use a similar approach to represent the soil described in chapter 2.

### 5.1.3 Unsaturated porous medium

To simplify the study, we assume a 1D porous medium, with a mean saturation in the porous layer of thickness $D$. We consider an unsaturated porous medium with a water saturation $s$ and a porosity $\phi$. We develop a gradient dynamics model to represent the water transfer as in the chapter 2. Neglecting the gravity free energy, the density of free energy of the system is the sum of two contributions: wetting energy including compression energy, and the free surface energy (see Eq.(2.25)). According to Eq.(2.43) the thermodynamic gradient dynamics model leads to

$$\phi \partial_t s = \nabla \cdot \left[ K(s) \nabla \left( \frac{\delta \mathcal{F}^\xi}{\delta s} - \rho g z \right) \right] \qquad (5.13a)$$

$$\mathcal{F}^\xi[s, c] = \int \left[ \frac{\Gamma}{2} |\nabla s|^2 + f_w^\xi(s, c) \right] dx^3, \qquad (5.13b)$$

$$\frac{\delta \mathcal{F}^\xi}{\delta s} = \partial_s f_w^\xi(s, c) - \Gamma \Delta s. \qquad (5.13c)$$

This model is similar to the droplet model on the previous Section 5.1.2. Note that the saturation $s$ and the height of the drop $\zeta$ have not the same dimension. Thus, we propose to define an equivalent height of water $\xi$ corresponding to this mean saturation $s$,

$$\xi = \phi D s. \qquad (5.14)$$
Therefore $\xi$ and $\zeta$ have the same dimension, i.e for uniform distributions of water height $\zeta$ and saturation, $\zeta + \xi$ represents the volume of water per unit area of our system that is constant throughout the dynamics. In other word,

$$ \int_{\Omega} (\zeta + \xi) d^2r = C, \qquad (5.15) $$

where $C$ is the mean height of water in the soil and on it.

Thus, inserting Eq.(5.14) in Eq.(5.13), we express the gradient dynamics model by

$$ \partial_t \xi = D \nabla \cdot \left[ K(\xi) \nabla \frac{\delta \mathcal{F}^\xi}{\delta \xi} \right], \qquad (5.16a) $$

$$ \mathcal{F}^\xi[\xi, c] = \int \left[ \frac{T}{2} |\nabla \xi|^2 + f_w^\xi(\xi, c) \right] dx^2, \qquad (5.16b) $$

$$ \frac{\delta \mathcal{F}^\xi}{\delta \xi} = \partial_\xi f_w^\xi(\xi, c) - T \Delta \xi. \qquad (5.16c) $$

where the surface tension coefficient $T = \frac{\Gamma}{\phi_D}$, the wetting energy in the soil $f_w^\xi(\xi, c) = \phi D f_w^\xi(s, c)$ and its derivative $\partial_\xi f_w^\xi(\xi, c) = -\rho g H_{cap} \Pi_w(\frac{\xi}{\phi_D}, c)$.

In the presence of high concentrations of amphiphilic matter $c$, the repellent energy is amplified. In other words, the repellent energy density $f_{rep}$ is related to an amplitude factor $\chi(c)$ depending on $c$ the concentration of amphiphilic matter. Thus, the dimensionless wetting pressure, $\Pi_w(s)$, associated to the wetting energy density $f_w^\xi$ is given by

$$ \Pi_w(s, c) = - \frac{\partial_s f_w^\xi}{\rho g H_{cap}} = h_c(s) + \chi(c) h_r(s) + h_{comp}(s), \qquad (5.17) $$

where $h_{comp}$ is the compressible pressure corresponding to the compressible energy. In this study, we adopt the Brooks and Corey function for the capillary pressure [42],

$$ h_c(s) = s^{-\frac{1}{\lambda}} \qquad (5.18) $$

and the compressible term $h_{comp}$ is

$$ h_{comp}(s) = \operatorname{sgn}(1 - \lambda) \exp(-\kappa(1 - s)) s^{-\frac{1}{\lambda}} \left(1 + \frac{\lambda}{\lambda - 1} \kappa s\right). \qquad (5.19) $$

The parameters $\lambda \neq 1$ and $\kappa \ge 0$ are respectively a pore size distribution and an air compression parameter [86]. Moreover, we fixed the repulsive effect to be present from $s_r = 0.05$ saturation and to be maximal around $s = 0.2$ with a decrease in its effects as saturation increases until it disappears at high saturation. These criteria lead us to adopt a Gaussian function centred at the repellent value $s_r$ with a dispersion coefficient $\sigma$,

$$ h_r(s) = - \frac{1}{\sigma \sqrt{2\pi}} \exp \left( - \frac{(s - s_r)^2}{2\sigma^2} \right). \qquad (5.20) $$

The resulting curve $\Pi_w(s, c)$ is shown in Figure 5.3. The function is strictly decreasing for low concentrations $c \le 0.1$. For higher concentration, the curve is no longer monotonic: $\Pi(s, c)$ increases
on the interval $[s_m, s_M]$. The saturation of $s_m$ and $s_M$ are respectively the saturation where repulsion is strongest and the saturation where hydrophobicity no longer acts. These values vary according to the value of $c$. The repellent effect growth at low saturation when the concentration of amphiphilic matter $c$ increase. For example, for $c = 0.8$ the disjoining pressure at $s_m(0.8) = 0.214$ is $\Pi_w(s_m(0.8), 0.8) = -2079$ Pa whereas for $c = 0$ at the same saturation the pressure is less influenced by repellent force, i.e the pressure is $\Pi_w(s_m(0.8), 0) = 13057$ Pa.

<figure><img src="image_88.png" /><figcaption>Figure 5.3: Wetting pressure in the soil for different amphiphilic concentration. Parameters are $\rho = 10^3$ kg.m$^{-3}$, $g = 9.81$ m.s$^{-2}$, $\mu = 10^{-3}$ Pa.s, $H_{cap} = 0.26$ m, $s_r = 0.05$, $\sigma = 0.2$, $\lambda = 0.972$, $\kappa = 200$, $\vartheta = 1.418$ m, $v = 1.402$ and $c_0 = 0.4$.</figcaption></figure>

For study the behavior of the porous medium with different concentration, we need to calibrate the amplitude factor $\chi(c)$. Indeed this factor represents a "repellent barrier" for a concentration greater than a threshold concentration [198]. In order to define it, we focus on the percolation experiment of Kroener based on capillary rise of water in stratified sand [198]. We chose a power law using an amplitude $\vartheta \ge 0$ and a constant exponent $v \ge 0$ to represent this amplitude factor,

$$ \chi(c) = \frac{\vartheta}{H_{cap}} c^{v}. \qquad (5.21) $$

As this calibration does not require the same dimension than $\zeta$, we used the system of equations (5.13) to represent the capillary rise through two layers of sand: one without amphiphilic matter and another with a concentration $c$ as described in Section 5.2.

### 5.1.4 Coupled drop/soil system

In this section, we study the coupling of the previous models if a drop is deposited on a soil.

#### 5.1.4.1 Interface Soil/Drop

Experiments in literature show that the spread of the drops depends on saturation [105, 219]. The drop spreading is wider on wet than on dry soil. In our model, this saturation-dependency of wettability implies that the free energy density $f^\zeta$ depends on $\xi$. Then, we assume that the repellent part of the
wetting drop density of free energy decreases as saturation increases [21]

$$f_w^\zeta(\zeta, \xi, c) = f_{wa}^\zeta(\zeta) + a(\xi, c)f_{wr}^\zeta(\zeta), \qquad (5.22)$$

where $f_{wa}^\zeta$ and $f_{wr}^\zeta$ are respectively the attractive and repellent wetting energies of the drop on a soil without amphiphilic matter. The amplitude $a$ depends on the water content, i.e. $\xi$ and on the concentration of amphiphilic matter $c$. According to the experiment, $a(\xi, c)$ is a decreasing function of $\xi$. We assume that the function $a(\xi, c)$ is written as

$$a(\xi, c) = p(c) \left( 1 - \frac{\xi}{\phi D} \right). \qquad (5.23)$$

The function $p(c)$ is an increasing function and is expressed as a parabolic function,

$$p(c) = c_0 + (1 - c_0)c^2. \qquad (5.24)$$

This expression of $p(c)$ reproduces the angle dependence observed empirically for a concentration $c$. The parameter $c_0$ is associated to the contact angle of a drop on an impermeable soil without amphiphilic matter.

Therefore, the wetting energy $f_w$ of the coupled system is given by

$$f_w(\zeta, \xi, c) = f_{wa}^\zeta(\zeta) + a(\xi, c)f_{wr}^\zeta(\zeta) + f_w^\xi(\xi, c). \qquad (5.25)$$

According to the general decomposition of the free energy Eq. (5.2), we identify the interface wetting energy $f^{int}$:

$$f_w^{int}(\zeta, \xi) = (a(\xi, c) - 1) f_{wr}^{\zeta}(\zeta), \qquad (5.26)$$

where $f_{wr}^\zeta$ is the repellent part of the drop wetting energy given in (5.11c).

This interface energy modifies the pressure in the drop but also in the soil by adding $\frac{\delta F^{int}}{\delta \zeta}$ and $\frac{\delta F^{int}}{\delta \xi}$, respectively:

$$\frac{\delta F^{int}}{\delta \zeta} = [a(\xi, c) - 1] \frac{\delta f_{wr}^{\zeta}}{\delta \zeta}(\zeta), \qquad (5.27)$$

$$\frac{\delta F^{int}}{\delta \xi} = \partial_{\xi} a(\xi, c) f_{wr}^{\zeta}(\zeta). \qquad (5.28)$$

The added pressure Eq. (5.27) is a hydrophilic interaction since it is always negative. This indicates that wet soil is more hydrophilic than dry soil. In a similar way, the added pressure Eq. (5.28) is always positive indicating that the soil is more hydrophobic when a thin water film is present. Such an interaction between a thin film and the soil is not described in the literature in the of soil framework. However, a similar interaction between drop and polymer brush is derived from a free energy approach in [323, 165, 155, 191].
### 5.1.4.2 Exchange of water

In this section, we model the water flux at the interface. We define the infiltration velocity $v$ of the drop to the porous medium proportional to the difference between the thin film and the soil pressure:

$$v \propto \left( \frac{\delta F}{\delta \zeta} - \frac{\delta F}{\delta \xi} \right)$$

As in the Washburn law in a porous medium [356, 142], we introduce a linear model to represent the infiltration flux:

$$v = \beta \left( \frac{\delta F}{\delta \zeta} - \frac{\delta F}{\delta \xi} \right) \qquad (5.29)$$

where $\beta > 0$ is called transfer rate.

This expression (5.29) of the infiltration velocity $v$ is consistent with dissipate free energy property (Appendix H.1).

Figure 5.4 shows the difference between the thin film and the soil pressure for uniform $(\zeta, \xi)$ profiles. It represents the water flux sign when only the wetting energy is considered. The flux of water goes from the film to the soil, if the flux density is positive. In this condition, that corresponds to the yellow area in Figure 5.4. In other words, the water in the drop phase infiltrates in the porous medium. For high saturation or an amphiphilic concentration above a critical concentration, we see the formation of a hydrophobic barrier (in blue in Figure 5.4). This formation corresponds to a region where flat film of water cannot infiltrate the medium. For example, for $c = 0.8$ a barrier is present for any film height when saturation in the soil is between $s = 0.155$ and $s = 0.32$.

### 5.1.4.3 Governing equations of the coupled drop/soil system

To summarize, for a drop infiltration on a porous medium, $\zeta$ and $\xi$ are governing by the coupled drop Eq.(5.12) and soil Eq.(5.16) including the interface terms Eq.(5.27) and Eq.(5.28):

$$\mathcal{F}[\zeta, \xi, c] = \int \left[ \frac{T}{2} |\nabla \xi|^2 + \frac{\gamma}{2} |\nabla \zeta|^2 + f_w^{\xi}(\xi, c) + f_w^{\zeta}(\zeta) + (a(\xi, c) - 1) f_{wr}^{\zeta}(\zeta) \right] dx^2, \qquad (5.30a)$$

$$\partial_t \zeta = \nabla \cdot \left[ Q(\zeta) \nabla \frac{\delta F}{\delta \zeta} \right] - \beta \left( \frac{\delta F}{\delta \zeta} - \frac{\delta F}{\delta \xi} \right), \qquad (5.30b)$$

$$\frac{\delta F}{\delta \zeta} = \partial_{\zeta} f_w^{\zeta}(\zeta) - \gamma \Delta \zeta + (a(\zeta, c) - 1) \partial_{\zeta} f_{wr}^{\zeta}(\zeta), \qquad (5.30c)$$

$$\partial_t \xi = D \nabla \cdot \left[ K(\xi) \nabla \frac{\delta F}{\delta \xi} \right] + \beta \left( \frac{\delta F}{\delta \zeta} - \frac{\delta F}{\delta \xi} \right), \qquad (5.30d)$$

$$\frac{\delta F}{\delta \xi} = \partial_{\xi} f_w^{\xi}(\xi, c) - T \Delta \xi + \partial_{\xi} a(\xi, c) f_{wr}^{\xi}(\xi). \qquad (5.30e)$$

To solve PDE system Eq.(5.30), we employ a finite element approach using the C++ FEM library OOMPH-LIB [169]. For time-stepping, a backward differentiation formula (BDF) scheme of second order is used, while both the time-steps and the spatial mesh are adaptive to the dynamics. The same method is used to solve Eq.(5.13).
<figure><img src="image_89.png" /><figcaption>Figure 5.4: Competition of pressure. Parameters are $ρ = 10^3$ kg.m$^{-3}$, $g = 9.81$ m.s$^{-2}$, $μ = 10^{-3}$ Pa.s, $b = 0.3679$, $ζ_p = 0.02$ mm, $Π_0 = 9514$ Pa, $H_{cap} = 0.26$ m, $s_r = 0.05$, $σ = 0.2$, $λ = 0.972$, $κ = 200$, $ϑ = 1.418$ m, $ν = 1.402$ and $c_0 = 0.4$.</figcaption></figure>

Before to study the coupled system, we need first to calibrate the parameters in the soil model (section 5.2).

## 5.2 Calibration of soil model

In this section, we focus on the amphiphilic porous media defined by Eq.(5.13). We aim to calibrate the function $χ(c)$ define in Eq.(5.21), representing the influence of the repellent part of the amphiphilic matter in a porous medium. In this section, we will focus on exopolysaccharides (EPS) in the rhizosphere, in particular on percolation experiments in stratified media with different concentrations of EPS.

In Kroener's work [198, 201], capillary rise experiments were carried out on sand samples with EPS similar to the rhizosphere. The studies were conducted on three ranges of media with different sand grain sizes. Their sandy media are classified as fine (0.125 – 0.2 mm in diameter), medium (0.2 – 0.355 mm) and coarse (0.63 – 1.0 mm). The 5 cm high samples are made up of a layer of approximately 5 mm of a mixture of EPS and sand representing the rhizosphere surrounded by two layers of dry untreated sand representing the average soil, without EPS. The samples were placed on suction tables with a hydraulic charge of 2.5 cm in order to fix a constant pressure at their base.

This experiment shows that after five minutes of capillary rise, depending on the concentration of EPS, two behaviors exist:

i) the water rises in the medium until it reaches its capillary height,
ii) the water does not rise through the layer containing EPS.

The situation ii) occurs beyond a threshold concentration $c_{crit}$. This critical concentration depends on the size grain of sand and more precisely decreases with coarser grain. This experimental result can be explained by a smaller surface area of coarse grains than of fine grains. Wettability energy is therefore lower as grain size increases. These experimental results provide a basis for calibrating the function $\chi(c)$ defined in Eq. (5.21).

In the numerical simulation, only two layers are considered: one representing the base made up of untreated sand with a height of 3.8 cm and the other the layer of mixture with a height of 1.7 cm. The domain length is fixed to $L_x = 55$ mm and we apply Dirichlet boundary condition at $z = 0$ and homogeneous Neumann boundary condition at $z = L_x$. The wettability of the stratified medium $\Pi(s, c, z)$ is governed by the wettability in each stratum defined by Eq. (5.17) as defined in Beltrame and Cajot [20]. It can be generalised by the relation,

$$ \Pi_w(s, c, z) = h_c(s) + f_{het}(z)\chi(c)h_r(s) + h_{comp}(s). \quad (5.31) $$

where $f_{het}$ is a heterogeneous function representing (in dashed line) in Figure 5.6 the region of influence of the repulsive part of the EPS,

$$ f_{het}(z) = \frac{1 + \tanh(100(z - 38.5))}{2}. \quad (5.32) $$

Here, the repellent interaction acts when $f_{het}(z) = 1$, i.e between 39 and 55 mm. Using Eq. (5.13) with the heterogeneous wettability in Eq. (5.31) and the data for the different media given in Kroener works [201, 199], we model the different range of media studied. The function $\chi(c)$ is presented in Figure 5.5. The experimental data from Kroener et al [198] are compared with our model as a coefficient influencing the repulsive forces. These data (red dot in Figure 5.5) were fitted using Eq. (5.21) (in blue line in Figure 5.5). As a result, the best fit gives us parameters $\vartheta = 1.418$ and $v = 1.402$.

<figure><img src="image_90.png" /><figcaption>Figure 5.5: Dependence of the repellent factor $\chi(c)$ on c.</figcaption></figure>

Numerical simulations with this fitted $\chi$ function reproduce capillary rise dynamics in different sandy media in the presence of different EPS concentrations, as shown in Figure 5.6. In fine sand, where the capillary height is $H_{cap} = 260$ mm and the pore size distribution $\lambda = 0.972$, the water rises in both strata for concentration for concentration lower than 1. The situation ii) occurs when the concentration is greater or equal to $c = 1$, as shown in Figure 5.6a. For a medium sand (Figure 5.6b) a
similar behavior occurs but for a lower value $c_{crit} = 0.75$ as predicted. For coarse grain, the situation is slightly different because the capillary height is $H_{cap} \approx 10$ mm is lower than the first layer height of sand. In Figure 5.6c, we see that the soil is almost saturated until the capillary rise then decreases gradually. At the interface mixture soil and soil without EPS, water is remaining at the saturation about $s = 0.5$. Above this interface, water further rises the EPS layer if concentration $c$ is below $c_{crit} = 0.05$.

<figure><img src="image_91.png" /><figcaption>Figure 5.6: Equilibrium profiles of saturation after capillary rise in sandy stratified media for three different grain sizes (a) fine grain, (b) medium grain and (c) coarse grain. The first layer ($0 \le z \le 39$ mm) does not contain EPS and the second layer ($z > 39$ mm) is a mixture of sand with EPS of different concentrations $c$ indicated in legend. The vertical dashed line delimits the two layers. The duration of capillary rise is indicated in the legend of each panel when water has infiltrated the second layer. Other parameters are $\rho = 10^3$ kg.m$^{-3}$, $g = 10$ m.s$^{-2}$, $s_r = 0.05$, $\sigma = 0.2$, $\kappa = 200$.</figcaption></figure>

Therefore, the model Eq. (5.13) of amphiphilic porous media is able to reproduce both regimes of capillary rise and the same dependence of $c_{crit}$ on grain size.

In the rest of this paper, we use this calibration also in the coupled model Eq. (5.30).

## 5.3 Equilibria of flat film and uniform saturation

In this section, we study the different equilibria of the system (5.30) on the particular case of flat film and uniform saturation. We aim to determine the stability of steady-states according to a water height $H$ such that

$$H = \xi + \zeta \qquad (5.33)$$

where $\zeta$ is the flat film height of water on the soil and $\xi$ the uniform height of water in the soil.

Note that this situation correspond to a coexistence of the flat film height and a uniform height of water in the soil. This description is similar than the sorption isotherm of a brush studied in Hartmann et al. [166].

As the closed thermodynamics system is dissipative, equilibria of the system are minima of the free energy with the constraint Eq.(5.33). The free energy decreases until an equilibrium profile $U_0 = (\zeta_0, \xi_0)$ which satisfies Eq.(5.33). This equilibrium profile is reached at an extremum of the function
$$ \Psi(\zeta, \xi, c) = f(\zeta, \xi, c) + \Lambda(H - \xi - \zeta), \qquad (5.34) $$

where $\Lambda$ is the Lagrange multiplier.

For any uniform height of saturation $\xi$ and flat film $\zeta$ with respect to the constraint Eq. (5.33), the wetting energy density $f_w$ represents the function $\Psi$, then equilibria $U_0$ are extrema of the wetting energy density. In other words, these steady-states are given by the system

$$ \left\{ \begin{aligned} \frac{\partial \Psi}{\partial \zeta}(U_0) &= \frac{\partial f_w}{\partial \zeta}(U_0) - \Lambda = 0, \\ \frac{\partial \Psi}{\partial \xi}(U_0) &= \frac{\partial f_w}{\partial \xi}(U_0) - \Lambda = 0, \\ \frac{\partial \Psi}{\partial \Lambda}(U_0) &= H - \xi - \zeta = 0, \end{aligned} \right. \qquad (5.35) $$

Note that this condition is equivalent to a no-flux exchange condition $v = 0$, i.e

$$ \frac{\delta F}{\delta \zeta}(U_0) - \frac{\delta F}{\delta \xi}(U_0) = 0. \qquad (5.36) $$

### 5.3.1 Equilibria for fixed $H$

Because of mass conservation, $H$ is constant. In this case, the wetting free energy is only expressed by the variable $\xi$. Let us define $\varphi(\xi)$ such as:

$$ \varphi(\xi) = f_w(\zeta(\xi), \xi), $$

where $\zeta(\xi) = H - \xi$ is the water film height.

Stable steady-states correspond to minima of $\varphi$ or in equivalent manner, the second derivative of $\varphi$ verifies (H.3),

$$ \frac{d^2\varphi}{d\xi^2}(\xi) = \partial_{\xi\xi}f_w - 2\partial_{\xi\xi}f_w + \partial_{\xi\xi}f_w > 0 \qquad (5.37) $$

In Figure 5.7, we represent the analysis of equilibria for a fixed total water height equal to the capillary length $H = \ell_c = \sqrt{\frac{\gamma}{\rho g}} \approx 2.65$ mm and a fine sandy soil with different EPS concentration $c$.

In the presence of amphiphilic matter, the free energy density curve (displayed in Figure 5.7a) is modified as a function of the concentration $c$. The solution $U_0$ of Eq. (5.35) changes with respect to $c$ and the initial profile in the soil. There are three case of equilibrium:

(i) There is only one global minimum which is reached when the medium is totally saturated, i.e. the compressed saturation $s_{comp} = 0.96$.

(ii) There is a global minimum corresponding to situation i) and an inflection point in the dry zone around a saturation $s_I$.

(iii) There is a global minimum corresponding to situation i) and two local extrema: a minimum in the dry zone around a saturation $s_m$ and a maximum around a saturation $s_M$.

These solutions are in agreement with the no-flux exchange condition in Figure 5.7b.
<figure><img src="image_92.png" /><figcaption>Figure 5.7: (a) Free energy density and (b) flux density with respect to saturation $s$. Parameters are the same as Figure 5.4 with $H = l_c$ mm and different concentration $c = \{0.1, 0.6, 0.71, 0.8, 1\}$.</figcaption></figure>

We find numerically, a critical concentration zone around $c_{crit} = 0.71$ separating situations i) and iii). We focus the study on the case i) and iii).

For low concentrations $c < c_{crit}$, corresponding to situation i), the unique solution is a minimum of the free energy $F$. When the concentration is greater than the critical concentration $c > c_{crit}$, there are three possible equilibrium states as described in situation iii). In this case, we have two additional equilibria corresponding to local extrema. These steady-states varying according to $c$. For example, the saturation of the minimum change from $s_m(0.8) = 0.151$ to $s_m(1) = 0.101$ whereas the maximum change from $s_M(0.8) = 0.317$ to $s_M(1) = 0.385$ (see Figure 5.7).

Only the two minima steady-states are stable whereas the maximum is unstable. To conclude, there exist two kinds of infiltration: a complete infiltration associated to the global minimum and a partial infiltration due to the repellent barrier (local minimum).

### 5.3.2 Equilibrium with respect to $H$

We study steady-states for different concentration and different total amount of water $H$ in the system. If initially, the porous medium is moist, i.e. the initial saturation is larger than $s_M(c)$, the equilibrium does not depend on $c$. In other words, equilibria $U_0$ after an infiltration into a wetted porous medium do not depend on the amphiphilic concentration $c$. These steady-states correspond to a minimum of the free energy density. Note that the profile $h_0$ increases linearly with $H$ once the saturation $s_0$ has reached a minimum energy.

In the case of dry initial porous medium, steady-states depend on the amphiphilic concentration $c$. Figure 5.8 displayed the simulated equilibrium with respect to $H$ and $c$.

For low height of water $H < \phi D$, there is not enough water to saturate the soil. As long as $c < 0.7$, the equilibrium saturation increases linearly with the height $H$ (Figure 5.8b).

Further increasing, the height of water $H > \phi D$, the soil is always saturated, for low concentration $c \le 0.71$. The energy density of these porous medium have only the global minimum as equilibrium. Since the soil is saturated, the height of the flat film increases linearly with the height $H$ as shown in Figure 5.8a.

Equilibrium saturation corresponds to the lower boundary of the hydrophobic barrier shown in blue in Figure 5.4d. This boundary is the local minimum defined by $s_m(c)$. Moreover, equilibrium height of
<figure><img src="image_93.png" /></figure>

<figure><img src="image_94.png" /><figcaption>Figure 5.8: Equilibria ($\zeta_0, s_0$) with respect to H for different concentration $c$.</figcaption></figure>

the flat film increases linearly with respect to the height $H$.

As for the fixed height of water, the model reproduces two infiltration kinds of uniform saturation and a flat film of water depending on the amphiphilic concentration. Stable equilibria are used to reproduce the initial soil profile before the drop is deposited. In the following, we focus on the infiltration of a drop in a soil in presence of EPS.

## 5.4 Drop infiltration

In this section, we numerically simulate the spreading and the infiltration of a water drop on a fine sandy soil (grain size about 200 µm) with an EPS concentration $c$. The volume per length of water drop is fixed $20.9 \text{ mm}^2$ and only the initial saturation of sandy soil and the EPS concentration $c$ vary.

### 5.4.1 Wet porous medium

We consider a pre-wetted fine sandy soil with a mean saturation $s = 0.8$ with a EPS concentration $c = 0.8$. According to the previous study the equilibrium state corresponds to an uniform saturation. We then place a drop on the surface ($20.9 \text{ mm}^2$ volume per length) and study its dynamics. Note that it can eventually infiltrate completely without saturating the soil.

The drop spreads rapidly over the substrate until it reaches the domain boundaries. This very fast dynamic takes about 0.1 s and the infiltration is negligible (see Fig. 5.9). The profile of the thin film of liquid is still drop-shaped: the drop height goes from 3 to 1 mm. Because of the boundary conditions, the water no longer spreads out, but the profile gradually flattens out over the course of a second. At the same time, the water gradually percolates down. Then the saturation reaches the capillary compressibility barrier at $t$ about 3 s. The dynamics converges to the same equilibrium computed in the previous section for the same total water height $H = 2.246 \text{ mm}$.

A wet sandy soil without EPS ($c = 0$) and the same initial saturation behaves in a similar way, both qualitatively and quantitatively. As a result, soil moisture hides the effects of the amphiphilic matter, and the soil surface behaves like a substrate with complete wetting properties.
<figure><img src="image_95.png" /><figcaption>(a) Drop phase.</figcaption></figure>

<figure><img src="image_96.png" /><figcaption>(b) Soil phase.</figcaption></figure>

Figure 5.9: Infiltration dynamic of a water drop, with a volume per length $V = 20.9 \text{ mm}^2$, in a wet soil, $s = 0.8$, and an EPS concentration $c = 0.8$. Parameters are the same as Figure 5.4.

### 5.4.2 Dry porous medium

We consider a dry fine sandy soil with different EPS concentrations $c$. We fix the saturation in the porous medium at $s = 0.1$. This value corresponds to a dry medium for which the uniform profile is in a stable zone for $c \le 1$ (see Figure 5.8).

<figure><img src="image_97.png" /><figcaption>(a) Drop phase.</figcaption></figure>

<figure><img src="image_98.png" /><figcaption>(b) Soil phase.</figcaption></figure>

Figure 5.10: Equilibria profiles of a drop on a dry sandy soil. Parameters are the same as Figure 5.4.

Equilibrium profiles of the liquid in each phase are displayed in Figure 5.10. These equilibria depends on concentration $c$.

For $c \le 0.7$ the drop infiltrates almost completely in the porous medium (see Figure 5.10a). The remaining volume fraction increases with concentration but remains below 12.6%. At low saturation ($c = 0.1$), the water infiltrates uniformly, whereas at higher saturation ($c \ge 0.6$), infiltration occurs essentially at the vertical of the drop. In the latter case, the equilibrium state is formed by a dry/wet structure corresponding to pseudo-partial diffusion (see Figure 5.10b). This type of structure was shown for hydrophobic soil in [20]. At $c = 0.8$ and larger concentration, the drop fails to infiltrate. It is stopped by the hydrophobic barrier formed in the soil. For instance, the saturation at equilibrium for a concentration of $c = 0.8$ is about $s = 0.1$ corresponding to the hydrophobic barrier boundary in Figure 5.4d.

Figure 5.11 shows that the infiltration rate increases with the concentration during the early stages. The imbibition duration increases with the concentration $c$. Moreover, depending on $c$ there are several stages associated to different infiltration rates. A break in slope occurs at $t = 1$ s for $c = 0.6$ resulting
in lower infiltration until total imbibition. At $c = 0.7$, these two stages occur again but the infiltration rate during the second stage is very low (between $t = 2$ and 14 s). A third stage (between $t = 14$ and 23 s) follows with a higher infiltration rate once again. Finally, a last, slower stage of infiltration occurs until total imbibition. Above a critical value of between $c = 0.7$ and $c = 0.8$, infiltration stops after a short time ($< 1$ s) and a drop setting on a porous medium remains, that defines an equilibrium.

<figure><img src="image_99.png" /><figcaption>Figure 5.11: Drop volume per length during infiltration. $V_{infiltration}$ in black dashed line corresponds to 12.6% of the initial drop that imbibed in the sandy soil. Parameters are the same as Figure 5.4.</figcaption></figure>

We detail in Figure 5.12 the drop and saturation profiles for $c = 0.7$ at different times stages. Initially, the drop spreads rapidly, $t = 10$ ms, over the substrate almost without infiltrating. Then, during the second stage, the radius of the drop decreases as it infiltrates vertically into the soil with little lateral diffusion. Infiltration is almost stopped as shown by the plateau of the infiltration curve in Figure 5.11. This behavior is due to the hydrophobic barrier in the hydrophobic saturation range around $s = 0.12$ as shown in Figure 5.12b. Once the hydrophobic saturation range has passed, infiltration is faster during the third and last stages than the second. At the end of the last stage, the drop is imbibated and the equilibrium state is reached after 35 s. There is a factor 1000 between the characteristic time of infiltration and the spreading. In order to study in detail the infiltration dynamics, we analyse the spreading radius of the drop, the evolution of the contact angle and the infiltration time.

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

<figure><img src="image_101.png" /><figcaption>Figure 5.12: Infiltration dynamic of a water drop, with a volume per length $V = 20.9 \text{ mm}^2$, in a dry soil, $s = 0.1$, and an EPS concentration $c = 0.7$. Parameters are the same as Figure 5.4.</figcaption></figure>

In the experimental literature, three infiltration stages were identified [228, 96, 177]:

i) Increasing Drawing Area (IDA): the radius of the drop increases while the angle decreases,
ii) Decreasing Drawing Area (DDA): the drop radius decreases and the contact angle remains constant,

iii) Constant Drawing Area (CDA): the drop radius remains constant due to the pinning of the contact line and the contact angle decreases.

Figure 5.13 represents the time evolution of drop radius and drop contact angle with logarithm scale to highlight the spreading dynamics. The numerical simulation of a soil with $c = 0.7$ displays the three stages dynamics described by i), ii) and iii) whereas for a soil with $c = 0.1$ and $c = 0.6$ display only the IDA and DDA dynamics (see Figure 5.13).

<figure><img src="image_102.png" /><figcaption>Figure 5.13: Infiltration step of a drop on a sandy soil. [Line] IDA, [Dashed] DDA and [dotted] CDA. Parameters are the same as Figure 5.4.</figcaption></figure>

In a sandy soil with $c = 0.1$, we have an IDA stage corresponding to a spreading behavior during 0.14 second. During this step, the contact angle goes from 60° to 11° and the drop radius from 4.19 to 11 mm. After the spreading stage, the drop starts a DDA stage where the infiltration drives the dynamic until the total imbibition in the substrate after 6.6 seconds. During this step, the contact angle decreases slightly until the equilibrium state corresponding to the complete imbibition.

From this study, the spreading times of the drop are lower than a millisecond for $c = 0.6$ and $c = 0.7$. It is noteworthy that these infiltration times are shorter than for soils with a concentration of $c = 0.1$. We believe that the hydrophobicity of the soil is the cause of this phenomenon, since it reduces drop spreading. In contrast, the infiltration time into the substrate is shorter for low concentrations, going from 6.6 seconds for $c = 0.1$ to 35 seconds for $c = 0.7$.

These trends are consistent with the experimentally observed dynamics. Our model is able to reproduce different experimental observation of infiltration dynamics in dry soil. We now examine infiltration on pre-wetted soil.

### 5.4.3 Pre-wetted porous medium

In this section, we simulate the experimental work that we have carried out and those of Charles-Williams [73]: deposit a second drop on a medium previously wetted by a first drop. Thus, we consider a soil previously wet by a first drop of 20.9 mm². The initial saturation profile in soil corresponds to the equilibrium profile found in Section 5.4.2 after imbibition of the first drop on dry soil.

The infiltration behavior of the second drop is similar to the first: if the drop has not infiltrated, the second drop does not infiltrate, otherwise it infiltrates. Thus, for $c > 0.8$, the drop infiltrates
slightly (Figure 5.14) whereas for $c = 0.7$, it infiltrates almost completely. At the equilibrium, we distinguish two saturation profiles in the soil. For low concentration $c = 0.1$, the saturation profile is almost uniform whereas for higher concentration $c = 0.6$ and $c = 0.7$, the water is confined in the range $x \in [-10 \text{ mm}, 10 \text{ mm}]$. This equilibrium is qualitatively the same for the first drop deposit (Section 5.4.2), only the wetted range in the soil is wider (see Figure 5.10). A similar behavior is pointed out in the chapter 2 and in Beltrame and Cajot [20] when the mean water content is increased.

<figure><img src="image_103.png" /><figcaption>Figure 5.14: Equilibria profiles of a drop on a pre-wetted sandy soil. Parameters are the same as Figure 5.4.</figcaption></figure>

In contrast the infiltration dynamics differs between the first and the second deposit of the drop. Only two stages of infiltration remain: first the IDA stage then the DDA stage. During the infiltration of the second drop, we have not the CDA stage for $c = 0.7$. This phenomenon is observed experimentally during the deposit of the second drop by Charles-Williams *et al.* [73]. As expected, the contact angle of the drop during infiltration is lower for pre-wetted soil than for dry soil (see Figure 5.15a). Comparing the time evolution of the drop volume between the first and second deposit, the infiltration rate is faster for the second deposit (see Figure 5.15b). However, after the infiltration plateau of the first

<figure><img src="image_104.png" /><figcaption>Figure 5.15: Infiltration dynamics of a drop on a dry and a pre-wetted sandy soil with $c = 0.7$. Parameters are the same as Figure 5.4.</figcaption></figure>

deposit, the infiltration rate is higher than for the second deposit. Thus the duration of imbibition is longer for the second deposit ($t = 50$ s) than for the first deposit on the dry substrate ($t = 35$ s).

Our model reproduces the longer times as well as the infiltration dynamics of the drop in a qualitative way between a dry soil and a wet soil.