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

HAL
open science

# Modeling water transfer in soil in the presence of amphiphilic matter : application to the rhizosphere

Florian Cajot

► To cite this version:

Florian Cajot. Modeling water transfer in soil in the presence of amphiphilic matter : application to the rhizosphere. Biological Physics [physics.bio-ph]. Université d'Avignon, 2024. English. NNT : 2024AVIG0623 . tel-04623150

HAL Id: tel-04623150

https://theses.hal.science/tel-04623150

Submitted on 25 Jun 2024

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

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

THÈSE DE DOCTORAT D'AVIGNON UNIVERSITÉ

École Doctorale n°536
Agrosciences et Sciences

Mention de doctorat :
Sciences Physiques

Laboratoire Environnement Méditerranéen & Modélisation des
Agro-Hydrosystèmes (EMMAH - UMR 1114)

Présentée par
Florian CAJOT

Modeling water transfer in soil in the presence of
amphiphilic matter: application to the rhizosphere

Soutenue publiquement le 1 mars 2024 devant le jury composé de :

<table><tr><td>Philippe BELTRAME</td><td>Maître de Conférences HDR</td><td>AU</td><td>Directeur de thèse</td></tr><tr><td>Annette BERARD</td><td>IPEF en chef HDR</td><td>INRAE</td><td>Co-encadrant</td></tr><tr><td>Claude DOUSSAN</td><td>Chargé de recherche</td><td>INRAE</td><td>Co-encadrant</td></tr><tr><td>Andrea SCHNEPF</td><td>Professeure</td><td>IBG Jülich</td><td>Rapporteuse</td></tr><tr><td>Uwe THIELE</td><td>Professeur</td><td>WWU Münster</td><td>Rapporteur</td></tr><tr><td>Laurent LIMAT</td><td>Directeur de Recherche</td><td>CNRS</td><td>Examinateur/Président du jury</td></tr><tr><td>Frédéric REES</td><td>Chargé de Recherche</td><td>INRAE</td><td>Examinateur</td></tr><tr><td>Cyprien SOULAINE</td><td>Chargé de Recherche HDR</td><td>ISTO</td><td>Examinateur</td></tr></table>
<figure><img src="image_3.png" /></figure>
Modeling water transfer in soil in the presence of amphiphilic matter: application to the rhizosphere
© 2024 by Florian CAJOT is licensed under Creative Commons Attribution 4.0 International. To
view a copy of this license, visit https://creativecommons.org/licenses/by/4.0/

<figure><img src="image_4.png" /></figure>
Aucune investigation humaine ne peux s'intituler véritable science,
s'il ne passe pas par la démonstration mathématiques

Léonardo da Vinci
Remerciements

Je tiens tout d'abord à remercier mon directeur de thèse Philippe Beltrame, Maître de Conférences à l'Université d'Avignon et mes co-encadrants Annette Bérard, IPEF en chef et Claude Doussan, Chargé de recherche à l'INRAE. Merci pour l'obtention de mon stage de M2, au sein de l'Unité Mixte de Recherche Environnement Méditerranéen et Modélisation des Agro-Hydrosystèmes (UMR, EMMAH), sans lequel je n'aurai pas pu réaliser ce travail de thèse. Je les remercie d'avoir suivi ces recherches, d'avoir partagé une part de leurs savoirs avec moi et de m'avoir conseillé tout au long de ces trois années et demi.

J'adresse mes sincères remerciements à l'INRAE et l'Université d'Avignon pour avoir financé ces travaux de recherche. Je profite également de l'occasion pour remercier l'Ecoles doctorales d'Avignon (ED 536) de m'avoir donné l'opportunité d'obtenir la bourse Perdiguierqui m'a permis d'établir et de maintenir des collaborations fructueuses avec les chercheurs de l'université de Munster en Allemagne.

J'adresse aussi toute ma gratitude à Laurent Limat, directeur de recherche, André Schnepf, Professeure à l'Université de Julich, Uwe Thiele, Professeur à l'Université de Munster, Frédéric Rees, chargé de recherche à l'INRAE et Cyprien Soulaine, chargé de recherche à l'ISTO qui m'ont fait l'honneur de faire partie du jury de thèse.

J'adresse également ma gratitude à Laurent Legendre, Laurent Limat, Stéphane Ruy, François Lafolie et Céline Pélosi pour leurs conseils et retour sur mes avancées scientifiques et mon parcours durant mes différents comités de pilotage / comités de suivi de thèse / comités de suivi individuel.

Merci aux collègues de Munster : Uwe, Simon, Tobias, Max et Fenna, vous avez été un tournant dans ma thèse. Vous m'avez permis de développer mes compétences en développement numérique et sur les concepts de thermodynamique.

Je remercie les collègues universitaires du département de physique pour leur aide en enseignement de physique et de m'avoir aider à transmettre mes connaissances aux étudiants ! Je remercie également les collègues de l'INRAE pour avoir répondu à toutes mes questions en physique du sol, de la faune et de la flore et des modélisations du sol ! Je remercie plus particulièrement Nicolas de m'avoir soutenu aussi bien dans les relectures de mon manuscrit, dans les rushs et pour nos nombreuses discussions avant la "fermeture du bâtiment".
Lucas, Marianne, Léa, Fédérico, Quentin et Léna, nos échanges ont enrichi mon expérience pluridisciplinaire ! Vous êtes les prochains sur la liste, courage à vous tous !

Clémentine, Emy et Achraf votre collaboration a été fondamentale pour l'élaboration d'une partie des résultats expérimentaux et numériques et pour cela je vous en remercie. Ce fut un plaisir de travailler avec vous !

Renaud, Laura, Harry, Uzair et Lucas, je vous remercie de votre bon humeur et des bons délires que nous avons partagé en tant que collègues de bureau. Je remercie particulièrement Léo pour tout le soutien que tu as pu m'apporter que ce soit durant les heures de boulot mais surtout en dehors. Nous avons partagé une belle aventure durant ton alternance et ma thèse !

Merci à Manon, Chloé, Hamza, Hamed, Elisabeth et Mostafa pour vos conseils durant cette épreuve de par votre expérience en tant que thésard et votre soutien que ce soit au travail et en dehors. Chloé, je te laisse le flambeau de la plus ancienne des non permanent de l'UMR en espérant de tout coeur que tu obtiennes un poste de titulaire !

Je souhaite remercier Clément et Cécile d'avoir été de super collègues et amis. Les nombreuses heures que nous avons passé ensemble à discuter de sciences et autres sujets du quotidien sont des moments chères à mes yeux. Votre aide et vos pensées critiques m'ont vraiment été d'un grand secours surtout avant mon rendu de manuscrit. Clément, courage pour le maraîchage et Cécile accroche toi tu vas réussir cette thèse !

Maxime, mon frère de coeur, ton soutien m'a été essentiel que ce soit durant cette thèse mais même bien avant ! Merci également pour ces moments en game avec Florentin ou nos parties de jeux de cartes qui ont été une véritable bouffée d'air frais durant cette épreuve !

Un grand remerciement également à mes amis de longue date : Chloé, Quentin, Florine et Baptiste pour tous les moments partagés ensemble !

Je voudrais aussi remercier ma famille : mes parents et mon frère. Vous m'avez soutenue durant cette thèse, comme vous l'avez fait tout au long de mes études.

Et pour finir, je remercie Marie ma partenaire d'aventure et binôme de vie. Malgré la distance et les rebondissements du quotidien, tu as été mon rayon de soleil et ma source de motivation, avec Lune. Merci pour ton amour, ton soutien et ta patience à toute épreuve. Je ne saurai t'exprimer toute ma gratitude alors je préfère te dédier ce manuscrit puisque tu en es le pilier principal.
Abstract

Organic matter (OM) is an important component of soils. In particular, in the rhizosphere (the soil volume influenced by root activity) specific OM compounds in the form of polymers such as ExoPolySaccharides (EPS) are present. These compounds, exuded by roots and micro-organisms, are amphiphile. They can influence the water transfer in the root zone and act as a buffer in root water uptake. Understanding this regulatory role of EPS and, more generally, how organic matter alters water transfer in the soil, are open questions in which modelling can give some insights.

To better account for the variable hydrophobic/hydrophilic nature of organic matter, we propose a paradigm shift in the representation of soil wettability. Instead of associating wettability with just a contact angle, we propose a conceptual energetic framework, in which an energy functional represents the effective interactions between the soil matrix and the free water surface. This energy of interaction contains both the attractive (hydrophilic) and repulsive (hydrophobic) interactions, in order to reproduce the amphiphilic nature of the molecules.

By including all the energies at work in the soil, i.e. energy associated with wettability, gravity and surface tension of the water, we can define the free energy of the system. We derived the governing equation of the water saturation leading to a Darcy-Richards type equation coupled with a free energy equation. This system of equations ensures consistency with thermodynamics.

The PDE system is solved numerically using finite element methods implemented in C++ with the Oomph-lib library based on adaptive spatial and temporal scheme. The code is applied to 1D, 2D and 3D axis-symmetries including different boundary conditions.

The model is first tested for its behaviour regarding hydrophobicity effects in soil. In this respect, it is able to reproduce the finite horizontal spreading of water over a hydrophobic porous medium. Applied to a vertically stratified porous medium with alternating hydrophilic and hydrophobic layers, the numerical simulation is able to reproduce, depending on the mean water content, the phenomena of water trapping in hydrophilic layers, the 'common' diffusive transfer (as the medium approaches saturation), intermittent flow and fingering instability. In addition, the numerical bifurcation analysis using the PDE2PATH library display hysteresis phenomena. This is the first time that a PDE system captures all of these phenomena in a unique soil model.

In the second part, we adapt and test from experiments our model to the specific characteristics of rhizospheric soil. We built a model medium consisting of sand with added EPS derived from chia mucilage (analogous of maize mucilage). Experiments were carried out to characterize the model parameters. In particular, we studied experimentally and numerically the infiltration of a drop placed on the sand/EPS artificial porous medium. Regarding the model, we show that the laws governing the water flow infiltrating the matrix are modified in the presence of organic matter. Once calibrated,
the model can reproduce various phenomena specific to the presence of amphiphilic matter observed
during our experiments or in the literature. In particular, the “buffering” effect during wetting and the
different stages of drop imbibition.

Finally, the model is applied to root water uptake, at the scale of a root segment with its rhizosphere,
exhibiting a concentration gradient of EPS with distance from root. This application enables us to
study the impact of the amphiphilic characteristics of the EPS in the rhizosphere during the water
extraction phase by the root or during the night rewetting phase of the rhizosphere, with or without
external water supply.

Integrated into a broader framework, at the scale of a root system, this modelling will contribute to
the analysis of rhizospheric traits likely to improve water use by plants in a context of climate change
and practices.

**Keywords of the thesis:** Porous media - Wettability - Exopolysaccharides - Free energy - Numerical modelling

<figure><img src="image_5.png" /></figure>
Résumé

La matière organique (MO) est un composant important des sols. En particulier, dans la rhizosphère (le volume de sol influencé par l'activité des racines), des composés spécifiques de la MO sont présents sous la forme de polymères tels que les ExoPolySaccharides (EPS). Ces composés, exsudés par les racines et les micro-organismes, sont amphiphiles. Ils peuvent influencer le transfert d'eau dans la zone racinaire et agir comme un tampon dans l'absorption de l'eau par les racines. Comprendre ce rôle régulateur des EPS et, plus généralement, comment la matière organique modifie le transfert d'eau dans le sol, sont des questions ouvertes pour lesquelles la modélisation peut apporter des éléments de réponse.

Pour mieux prendre en compte la nature variable hydrophobe/hydrophile de la matière organique, nous proposons un changement de paradigme dans la représentation de la mouillabilité du sol. Au lieu d'associer la mouillabilité à un simple angle de contact, nous proposons un cadre conceptuel énergétique, dans lequel une fonctionnelle énergétique représente les interactions effectives entre la matrice du sol et la surface de l'eau libre. Cette énergie d'interaction contient à la fois les interactions attractives (hydrophiles) et répulsives (hydrophobes), afin de reproduire la nature amphiphile des molécules.

En incluant toutes les énergies à l'œuvre dans le sol, c'est-à-dire l'énergie associée à la mouillabilité, à la gravité et à la tension superficielle de l'eau, nous pouvons définir l'énergie libre du système. Nous avons dérivé l'équation régissant la saturation en eau, ce qui conduit à une équation de type Darcy-Richards couplée à une équation d'énergie libre. Ce système d'équations assure la cohérence avec la thermodynamique.

Le système d'EDP est résolu numériquement à l'aide de méthodes d'éléments finis implémentées en C++ avec la bibliothèque Oomph-lib basée sur un schéma spatial et temporel adaptatif. Le code est appliqué à des systèmes 1D, 2D et 3D axisymétriques incluant différentes conditions aux limites.

Le modèle est d'abord testé pour son comportement concernant les effets d'hydrophobie dans le sol. À cet égard, il est capable de reproduire l'étalement horizontal fini de l'eau sur un milieu poreux hydrophobe. Appliquée à un milieu poreux stratifié verticalement avec une alternance de couches hydrophiles et hydrophobes, la simulation numérique est capable de reproduire, en fonction de la teneur moyenne en eau, les phénomènes de piégeage de l'eau dans les couches hydrophiles, le transfert diffusif "commun" (à mesure que le milieu approche de la saturation), l'écoulement intermittent et l'instabilité de doigté. En outre, l'analyse numérique de bifurcation utilisant la bibliothèque PDE2PATH met en évidence des phénomènes d'hystérésis. C'est la première fois qu'un système EDP capture tous ces phénomènes dans un modèle de sol unique.

Dans la deuxième partie, nous adaptons et testons expérimentalement notre modèle aux caractéristiques spécifiques du sol rhizosphérique. Nous avons construit un milieu modèle composé de sable
additionné d'EPS dérivé du mucilage de chia (analogue au mucilage de maïs). Des expériences ont été
réalisées pour caractériser les paramètres du modèle. En particulier, nous avons étudié expérimentale-
ment et numériquement l'infiltration d'une goutte placée sur le milieu modèle sable/EPS. Concernant
le modèle, nous montrons que les lois régissant le flux d'eau s'infiltrant dans la matrice sont modifiées
en présence de matière organique. Une fois calibré, le modèle peut reproduire différents phénomènes
spécifiques à la présence de matière amphiphile observés lors de nos expériences ou dans la littérature.
En particulier, l'effet "tampon" lors du mouillage et les différentes étapes de l'imbibition des gouttes.

Enfin, le modèle est appliqué à l'absorption d'eau par les racines, à l'échelle d'un segment de
racine avec sa rhizosphère, présentant un gradient de concentration en EPS en fonction de la distance
à la racine. Cette application nous permet d'étudier l'impact des caractéristiques amphiphiles des
EPS dans la rhizosphère pendant la phase d'extraction d'eau par la racine ou pendant la phase de
réhumidification nocturne de la rhizosphère, avec ou sans apport d'eau externe.

Intégrée dans un cadre plus large, à l'échelle d'un système racinaire, cette modélisation contribuera
à l'analyse des traits rhizosphériques susceptibles d'améliorer l'utilisation de l'eau par les plantes dans
un contexte de changement climatique et de pratiques.

Mots-clés de la thèse : Milieu poreux - Mouillabilité - Exopolysaccharides - Energie libre - Mod-
élisation numérique

<figure><img src="image_6.png" /></figure>
Nomenclature

• $V_t$: Total volume of porous medium [m$^3$]

• $V_f$: Fluid volume [m$^3$]

• $V_s$: Solid volume [m$^3$]

• $V_a$: Air volume [m$^3$]

• $V_w$: Water volume [m$^3$]

• $\phi$: Porosity [–]

• $I_v$: Void index [–]

• $V_b$: Specific volume [m$^3$.kg$^{-1}$]

• $M_d$: Dry mass [kg]

• $\rho_b$: Apparent density of a porous medium in a dry state [kg.m$^{-3}$]

• $w$: Gravimetric water content [kg.kg$^{-1}$]

• $M_w$: Mass of water [kg]

• $M_d$: Mass of dry solid particles [kg]

• $\theta$: Volumetric water content [m$^3$.m$^{-3}$]

• $s$: Saturation of water [m$^3$.m$^{-3}$]

• $\gamma_d$: Dry density [–]

• $\rho_d$: Density of the dry medium [kg.m$^{-3}$]

• $\rho_w$: Density of the water [kg.m$^{-3}$]

• $\theta_{sat}$: Saturated volumetric water content [m$^3$.m$^{-3}$]

• $\theta_e$: Effective volumetric water content [m$^3$.m$^{-3}$]

• $\theta_r$: Residual volumetric water content [m$^3$.m$^{-3}$]

• $s_e$: Effective saturation of water [m$^3$.m$^{-3}$]
• $s_a$: Saturation of air [m$^3$.m$^{-3}$]

• $s_r$: Saturation of residual water (bounded to the solid) [m$^3$.m$^{-3}$]

• $s_w$: Saturation of free water [m$^3$.m$^{-3}$]

• $s_{comp}$: Compression saturation threshold [m$^3$.m$^{-3}$]

• $s_m$: mean saturation [m$^3$.m$^{-3}$]

• $\Phi$: Hydraulic potential per unit mass [J.kg]

• $\Psi$: Hydraulic potential per unit volume of water [J]

• $\phi$: Hydraulic potential per unit weight of water [J.N$^{-1}$]

• $Re$: Reynolds number [–]

• $Q$: Volumetric flow rate [m$^3$.s$^{-1}$]

• $S$: Cross-section of a porous medium [m$^2$]

• $L$: Length of the porous medium [m]

• $\mu$: Dynamic viscosity of the fluid [Pa.s]

• $k$: Intrinsic permeability of a porous medium [m$^2$]

• $K_{sat}$: Saturated hydraulic conductivity [m.s$^{-1}$]

• $g$: Gravity [m.s$^{-2}$]

• $P$: Pressure [Pa]

• $H$: Pressure head [m]

• $K$: Hydraulic conductivity [m.s$^{-1}$]

• $\Phi_t$: Total hydraulic potential [J.kg$^{-1}$]

• $h(\theta)$: Water retention characteristic curve [m]

• $\rho_p$: Density of fluid in a porous medium [kg.m$^{-3}$]

• $C(\theta)$: Specific capillary capacity [m$^{-1}$]

• $D(\theta)$: Specific diffusion coefficient [m$^2$.s$^{-1}$]

• $\gamma$: Surface tension [N.m$^{-1}$]

• $W$: work [J]

• $A$: surface [m$^2$]

• $F$: Helmholtz free energy [J]
• $l_c$: capillary length [m]

• $S$: spreading parameter [J.m$^{-2}$]

• $\Theta$: contact angle [rad]

• $H_{cap}$: capillary length [m]

• $\zeta$: thickness film [m]

• $\Pi$: disjoining pressure [Pa]

• $\zeta_p$: precursor film [m]

• $J$: Landau free energy [J]

• $f$: density free energy (for saturation [J.m$^{-3}$] and for thin film [J.m$^{-2}$])

• $\Gamma$: macroscopic surface tension coefficient [N]

• $\mathcal{F}$: hydraulic flux [m.s$^{-1}$]

• $R$: Rayleighian

• $\frac{\delta \mathcal{F}}{\delta s}$: variational derivative of $\mathcal{F}$ [J.m$^{-3}$]

• $f_{het}$: heterogeneous function [–]

• $L^*$: dimensionless length [–]

• $T^*$: dimensionless time [–]

• $\chi$: repellent amplitude parameter

• $r_{root}$: root radius [m]

• $r_{rhizo}$: rhizosphere radius [m]

• $r_{bulk}$: bulk soil radius [m]
# Contents

Introduction 1

I State of art: The basics of soil water transfer models 3

1 The basics of soil water transfer models 5

1.1 Description of a soil 6

1.1.1 Characteristic parameters and notations of a porous medium 6

1.1.2 Organic matter: a substance influencing soil structure 10

1.1.3 The hydraulic properties of a soil: application to the rhizosphere 12

1.2 Model of water transfer 16

1.2.1 Model without amphiphilic matter 16

1.2.2 Model in effective medium in presence of mucilage 22

1.2.3 Pore scale model 25

1.3 Modeling of wettability and application to the soil 25

1.3.1 Wettability on plane substrate 25

1.3.2 Wettability in the soil 30

1.3.3 Wettability energies 32

1.4 Conclusion 35

II Model and numerical scheme 37

2 Wettability and free energy in soil 39

2.1 Free energy in soil 40

2.1.1 Thermodynamic approach 40

2.1.2 Free energy approach 42

2.1.3 Bulk energy 44

2.1.4 Wettability Energy 44

2.1.5 Water free surface energy 48

2.1.6 Conclusion 50

2.2 Water transfer model 50

2.2.1 Homogeneous porous medium 51

2.2.2 Stratified porous medium 53
<table><tbody><tr><td>2.2.3 Boundaries conditions</td><td>55</td></tr><tr><td>2.2.4 Dimensionless water transfer in homogeneous medium</td><td>56</td></tr><tr><td>2.3 Theoretical results of homogeneous soil</td><td>57</td></tr><tr><td>2.3.1 Wetting pressure</td><td>57</td></tr><tr><td>2.3.2 Equilibrium</td><td>58</td></tr><tr><td>2.3.3 Linear stability of the uniform steady-state</td><td>61</td></tr><tr><td>2.4 Instabilities in stratified media</td><td>64</td></tr><tr><td>2.4.1 Conclusion</td><td>65</td></tr><tr><td>2.5 Publication of the model</td><td>65</td></tr><tr><td>2.5.1 Letter journal: EPL</td><td>66</td></tr><tr><td>2.5.2 Summary of the article</td><td>66</td></tr><tr><td>2.5.3 Article</td><td>67</td></tr><tr><td>2.6 Conclusion</td><td>75</td></tr><tr><td><b>3 Development and use of numerical tools to study water transfers</b></td><td><b>77</b></td></tr><tr><td>3.1 State of the art of the time integration</td><td>78</td></tr><tr><td>3.1.1 Numerical methods</td><td>78</td></tr><tr><td>3.1.2 Numerical time schemes</td><td>79</td></tr><tr><td>3.1.3 Numerical problem statement</td><td>79</td></tr><tr><td>3.1.4 My approach</td><td>81</td></tr><tr><td>3.2 Finite element method</td><td>82</td></tr><tr><td>3.2.1 Weak formulation</td><td>82</td></tr><tr><td>3.2.2 Galerkin method</td><td>82</td></tr><tr><td>3.3 Meshing</td><td>84</td></tr><tr><td>3.3.1 Spatial mesh</td><td>84</td></tr><tr><td>3.3.2 Time mesh</td><td>87</td></tr><tr><td>3.4 Numerical temporal scheme</td><td>91</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>3.5 Numerical Bifurcation Analysis</td><td>95</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>3.6 Conclusion</td><td>100</td></tr><tr><td><b>III Application to the rhizosphere</b></td><td><b>101</b></td></tr><tr><td><b>4 Characterization and parameterization of soil water properties in the presence of amphiphilic matter</b></td><td><b>103</b></td></tr><tr><td>4.1 An artificial rhizosphere medium</td><td>104</td></tr><tr><td>4.1.1 Choosing of root mucilage for building rhizosphere soil</td><td>104</td></tr><tr><td>4.1.2 The porous media design for rhizosphere analogue</td><td>109</td></tr></tbody></table>
<table><tbody><tr><td>4.1.3</td><td>Design and elaboration of the analogue rhizosphere soil</td><td>110</td></tr><tr><td>4.2</td><td>Characterization of hydraulic properties of the rhizosphere analogue medium</td><td>112</td></tr><tr><td>4.2.1</td><td>Materials and methods</td><td>113</td></tr><tr><td>4.2.2</td><td>Results</td><td>118</td></tr><tr><td>4.3</td><td>Parametrization of soil hydraulic properties accounting for variable wettability</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><tr><td>5</td><td>Water transfer model in the rhizosphere in the presence of EPS</td><td>129</td></tr><tr><td>5.1</td><td>Model of drop infiltration in rhizosphere analogue soil</td><td>130</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.2</td><td>Calibration of soil model</td><td>138</td></tr><tr><td>5.3</td><td>Equilibria of flat film and uniform saturation</td><td>140</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><tr><td>5.4</td><td>Drop infiltration</td><td>143</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>5.5</td><td>Water transfer around a single root</td><td>148</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>5.6</td><td>Humectation of the soil</td><td>151</td></tr><tr><td>5.7</td><td>Conclusion</td><td>154</td></tr><tr><td colspan="2">Conclusion</td><td>155</td></tr><tr><td colspan="2">Appendices</td><td>161</td></tr><tr><td colspan="2">A Numerical time scheme</td><td>161</td></tr><tr><td>A.1</td><td>Time scheme</td><td>161</td></tr><tr><td>A.1.1</td><td>Time scheme at one step</td><td>162</td></tr><tr><td>A.1.2</td><td>Multi-step time scheme</td><td>163</td></tr><tr><td>A.2</td><td>Prediction-correction methods</td><td>165</td></tr><tr><td>A.2.1</td><td>General principle</td><td>165</td></tr><tr><td>A.2.2</td><td>Order 1</td><td>166</td></tr><tr><td>A.2.3</td><td>Order 2</td><td>166</td></tr></tbody></table>
<table><tbody><tr><td>B Matrix time system</td><td>169</td></tr><tr><td>B.1 Linear operator</td><td>169</td></tr><tr><td>B.2 Backward Euler scheme</td><td>170</td></tr><tr><td>B.3 Trapezoid rule scheme</td><td>170</td></tr><tr><td>B.4 BDF2 scheme</td><td>171</td></tr><tr><td>C Variational formulation of the problem</td><td>173</td></tr><tr><td>C.1 Boundary conditions</td><td>174</td></tr><tr><td>C.1.1 General case</td><td>174</td></tr><tr><td>C.1.2 In dimension 1</td><td>175</td></tr><tr><td>C.1.3 NH condition on free energy input, and a flux condition on output</td><td>177</td></tr><tr><td>D Numerical resolution in time</td><td>179</td></tr><tr><td>D.1 Picard method</td><td>180</td></tr><tr><td>D.2 Newton method</td><td>180</td></tr><tr><td>D.3 Synthesis</td><td>182</td></tr><tr><td>E Chemical diversity of crop plants root mucilages. Implications for their maximal water content and mineralization in soils</td><td>185</td></tr><tr><td>F Material and methods : Designing a rhizosphere soil</td><td>199</td></tr><tr><td>G Experiments on the degradation of samples of chia mucilage mixed with sand</td><td>205</td></tr><tr><td>G.1 Materials and methods:</td><td>205</td></tr><tr><td>G.1.1 Microbial degradation</td><td>205</td></tr><tr><td>G.1.2 Presence of ExoPolySaccharides (EPS) adhering to the sand</td><td>206</td></tr><tr><td>G.2 Results and discussion</td><td>206</td></tr><tr><td>H Coupled system: Drop infiltration in porous medium with amphiphilic matter</td><td>209</td></tr><tr><td>H.1 Dissipative free energy</td><td>209</td></tr><tr><td>H.2 Model parameters</td><td>210</td></tr><tr><td>H.3 Stability analysis</td><td>211</td></tr><tr><td>H.3.1 Uniform variation</td><td>212</td></tr><tr><td>H.3.2 Equilibrium stability</td><td>212</td></tr><tr><td>H.4 Extrema of free energy density</td><td>214</td></tr><tr><td>I Water transfer with single root</td><td>215</td></tr><tr><td>I.1 Finite Element Method</td><td>215</td></tr><tr><td>I.1.1 Dimensional formulation</td><td>215</td></tr><tr><td>I.1.2 Dimensionless formulation</td><td>216</td></tr><tr><td>I.2 Simulated parameters</td><td>216</td></tr><tr><td>J Publications and Seminars</td><td>219</td></tr><tr><td>J.1 Publications in peer-reviewed journals</td><td>219</td></tr><tr><td>J.2 Seminars and conferences</td><td>219</td></tr></tbody></table>
# Introduction

## General context

The effects of greenhouse gases in the atmosphere are rapidly warming the earth's surface [309]. This rise in average temperatures is having a major impact on oceans, glaciers, hydrological cycles and biodiversity [118].

Among the impacts of climate change, average changes in annual hydrology have been observed. In terms of precipitation, some zones are impacted by an increase in average precipitation, such as the equator and high latitudes, while others are suffering from heatwaves [309]. In Central Europe, we are seeing increasingly hot summers with severe heatwaves [10], particularly in the Mediterranean zone.

Rising temperatures and climatic deregulation over extended periods impact soil and plant cover (hydric stress influences the reproduction and performance of plants [181, 308]). This deregulation influences water resources in the soil, such as groundwater recharge. In dry periods, it is difficult for water to infiltrate into the soil, and the small amount of water that infiltrates is used to irrigate crops, which has a major impact on water resources [192].

Before trying to find ways of improving plants' resistance to water stress, we need to understand the mechanisms used by plants and soils to influence water transfer, such as the release of organic matter (mucilage for plants and Extracellular Polymeric Substance for microorganisms [185]). There is a lot to learn about the underground part of plants, which is a much less well known zone than the aerial part.

The release of these amphiphilic substances<fnref n="1" /> must have a role to play for soil organisms, since they require an energy cost from them that could have been allocated to other functions. Therefore, it is interesting to study the influence of these substances in the soil, and in particular on hydraulic transfer such as water uptake by roots in the rhizosphere. Almost all the water used by the plant passes through zones of the soil influenced by these substances. This water is uptaken by the plant's root system to satisfy the environmental demands of the aerial part (evapotranspiration process).

This zone of soil close to the roots is distinguished from bulk soil by its special relationship with the plant: it is called the rhizosphere [306]. However, this zone is complex to study because of its location in the soil. Mathematical models are therefore needed to study the processes and simulate the hydraulic behavior of the rhizosphere under the influence of amphiphilic matter.

<footnote n="1">¹chemical substances containing a part with an affinity for water (hydrophilic) and an another one which repels it (hydrophobic).</footnote>
# The objectives

The aim of this thesis is to go further than current approaches to model water transfer in soils by taking into account the wettability of organic matter in the soil.

In this work, we will pay particular attention to a thermodynamic concept affecting soil dynamics: the free energy of water in the soil. Thus, knowing the hydraulic properties of the soil makes it possible to couple processes in it, then better understand the processes associated with amphiphilic matter and water in the soil.

We developed a model for simulating water transfer in the rhizosphere. Due to the high non-linearity of this new model, we will propose tools for temporal analysis and simulation. The objectives of this thesis are:

* improve water transfer models in soils, particularly for macroscopic models taking into account the presence of amphiphilic matter;

* provide numerical simulation tools;

* understand the influence of organic matter in the soil (e.g. EPS in the rhizosphere).

# Outline

This manuscript is divided into three main parts.

The first part presents a state of the art on the issue of organic matter influence in soil and the modeling of their function in water transfer. The objectives of this chapter 1 are (i) to understand the composition of a soil and the effects of organic matter on its hydraulic properties, and (ii) to understand the modeling of water transfer in a soil in the presence of organic matter such as in the rhizosphere. Moreover, this chapter suggests ways of improving the modeling of water transfer taking into account the concept of wettability.

The second part is devoted to the modeling of water transfer in a porous medium in the presence of amphiphilic matter. The chapter 2 presents the analysis and theoretical results of water transfer modeling based on the concept of free energy of water in soil. We will see in this chapter (i) how to represent the free energy of water in a soil and (ii) how wettability in the soil affects water transfer dynamics. The chapter 3 concerns numerical tools for simulating water transfer. We will present the numerical difficulties in simulating and analyzing transfer dynamics with this new model.

The third and final part of the manuscript is devoted to the application of this modeling in the rhizosphere. The chapter 4 presents the experimental approaches for characterizing and parameterizing an analogue rhizosphere medium. This chapter shows (i) how to design an analogue rhizosphere medium and (ii) how to characterize and parameterize its hydraulic properties. Finally, the chapter 5 presents applications of water transfer dynamics in the rhizosphere, such as the infiltration of a drop into an analogue rhizosphere medium, or water uptake by a root.

Finally, this manuscript ends with a general conclusion and possible perspectives for extending this work.
# Part I

## State of art: The basics of soil water transfer models
<figure><img src="image_7.png" /></figure>
# Chapter 1

## The basics of soil water transfer models

The aim of this chapter is to define the theoretical concepts associated to water transfer in soils, and to present existing work in the literature on the amphiphilic matter that influences water transfer dynamics in soil. We review the description of a soil and the effects of organic matter in the soil (section 1.1), to describe water transfer models in the soil and more particularly in the rhizosphere (section 1.2). To complete the global view of this thesis, we present the notion of wettability on a substrate as a perspective for improving models of water transfer in section 1.3.

### Contents

<table><tr><td>1.1</td><td>Description of a soil</td><td>6</td></tr><tr><td>1.1.1</td><td>Characteristic parameters and notations of a porous medium</td><td>6</td></tr><tr><td>1.1.2</td><td>Organic matter: a substance influencing soil structure</td><td>10</td></tr><tr><td>1.1.3</td><td>The hydraulic properties of a soil: application to the rhizosphere</td><td>12</td></tr><tr><td>1.2</td><td>Model of water transfer</td><td>16</td></tr><tr><td>1.2.1</td><td>Model without amphiphilic matter</td><td>16</td></tr><tr><td>1.2.1.1</td><td>Saturated medium: Darcy law</td><td>16</td></tr><tr><td>1.2.1.2</td><td>Unsaturated porous medium</td><td>17</td></tr><tr><td>1.2.2</td><td>Model in effective medium in presence of mucilage</td><td>22</td></tr><tr><td>1.2.3</td><td>Pore scale model</td><td>25</td></tr><tr><td>1.3</td><td>Modeling of wettability and application to the soil</td><td>25</td></tr><tr><td>1.3.1</td><td>Wettability on plane substrate</td><td>25</td></tr><tr><td>1.3.2</td><td>Wettability in the soil</td><td>30</td></tr><tr><td>1.3.3</td><td>Wettability energies</td><td>32</td></tr><tr><td>1.4</td><td>Conclusion</td><td>35</td></tr></table>
# 1.1 Description of a soil

Soil is a complex porous medium containing minerals and organisms (living and dead). It is represented by three main phases:

* the solid phase composed of organic matter and minerals (sand, silt and clay);

* the liquid phase composed of inorganic ions and ionized and non-ionized organic molecules;

* the gaseous phase close to air at the soil surface, but which may be enriched at depth with gaseous compounds such as nitrogen oxides(NO,N₂O).

In this section, we present the complexity of a soil, describing the characteristics of a porous medium, such as the structure of a soil in the presence of organic matter, in order to define the hydraulic properties of a soil.

## 1.1.1 Characteristic parameters and notations of a porous medium

A porous medium has physical characteristics associated with dynamics such as the flow of water through it. Characteristic parameters are therefore needed to describe and model the structure and porosity of the medium, or its moisture content. These characteristic parameters are essential for representing the properties of the medium in models and simulations.

### Structure and porosity

A porous medium is formed by solid elements and void space (Figure 1.1). Here we refer to the solid phase as the solid matrix, and the void phase as the pore space or fluid space. The total volume of the porous medium $V_t$ [m³] is the sum of fluid volume $V_f$ [m³] and solid volume $V_s$ [m³],

$$V_t = V_f + V_s. \qquad (1.1)$$

The pore space can contain one or more fluids (e.g. water and air). In the case of water transfer in

<figure><img src="image_8.png" /><figcaption>Figure 1.1: Sketch of a porous medium, inspired of [91].</figcaption></figure>

the medium, the pore space of volume $V_f$ comprises a volume of air $V_a$ [m³] and a volume of water $V_w$ [m³] such that

$$V_f = V_a + V_w. \qquad (1.2)$$
The pore space depends on its available fluid volume. Indices of void volume relative to the structure of the medium are used to characterize the fluid space. Among these indices, we have the following:

* porosity $\phi$ [–] is a soil physics notation representing the index of relative pore void volume $V_f$ in soil $V_t$,

$$ \phi = \frac{V_f}{V_t}; \qquad (1.3) $$

* void index $I_v$ [–] is a soil mechanics and engineering notation representing the index of relative void volume of pores $V_f$ compared to solid volume $V_s$,

$$ I_v = \frac{V_f}{V_s}. \qquad (1.4) $$

In the rest of this manuscript, we will use the porosity notation to represent the fluid space. Note that this $\phi$ parameter is between 0 and 1.

Porosity size varies in porous media. A classification of porosity in the medium is based on the void size of the pore space. This criterion's limitations in classifying porosity are much debated in the literature [14, 31]. Luxmoore [222] suggests a classification that associates a type of porosity with a void size:

* macroporosity: size greater than 100 µm;

* mesoporosity: size between 10 and 100 µm;

* microporosity: size smaller than 10 µm.

However, porosity also depends on the solid structure of the medium. The distribution of solid elements influences the pore space. The more compact the medium, the smaller the void volume is. We then express a specific volume $V_b$ [m³.kg⁻¹] indicating the compactness index representing the total volume of the medium in relation to its dry mass $M_d$ [kg],

$$ V_b = \frac{V_s}{M_d}, \qquad (1.5) $$

where $M_d$ is the mass of dry solid particles [kg]. Note that this specific volume is the inverse of the apparent density of a porous medium in a dry state $\rho_b$ [kg.m⁻³],

$$ \rho_b = \frac{M_d}{V_t}. \qquad (1.6) $$

Apparent density and specific volume are two indicators of the compaction of a porous medium. The bulk density is the most widely used compaction indicator in the literature.

## Moisture of the porous medium

Water is a fluid commonly found in porous media. This liquid flows into the pore space. A level of moisture is then present in porous medium. There are several indicators that express this moisture:
* weight moisture or gravimetric water content $w$ [kg.kg⁻¹] is the mass of water $M_w$ [kg] relative to the mass of dry solid particles $M_d$,

$$w = \frac{M_w}{M_d}; \qquad (1.7)$$

* volumetric water content $\theta$ [m³.m⁻³] is the volume of water $V_w$ in the volume of the medium $V_t$,

$$\theta = \frac{V_w}{V_t}; \qquad (1.8)$$

* degree of saturation or saturation $s$ [m³.m⁻³] is the volume of water $V_w$ in the volume of pores $V_f$,

$$s = \frac{V_w}{V_f}. \qquad (1.9)$$

Note that these moisture indicators are equivalent. Each indicator is appropriate to the moisture study to be carried out. Gravimetric water content is very useful for measuring moisture in porous media. However, it is not representative of moisture in terms of the volume of water available and present in the medium. In this situation, volumic water content and saturation indicators are more appropriate. These indicators are used to study the moisture in the volumes of the medium. There exist equivalence relationships between each moisture. Among these relationships, we have these main equivalences [128, 170]:

* the gravimetric water content $w$ and the volumic water content $\theta$ are related by the dry density $\gamma_d$ [-] of the medium and is written as

$$\theta = w\gamma_d, \qquad (1.10)$$

where $\gamma_d$ is the ratio of the density of the dry medium $\rho_d$ [kg.m⁻³] to the density of the water $\rho_w$ [kg.m⁻³],

$$\gamma_d = \frac{\rho_d}{\rho_w}; \qquad (1.11)$$

* saturation $s$ is expressed by the ratio between the volume weight of the gravimetric water content $\rho_dgw$ [kg.m⁻³.kg⁻¹.m³] and the volume weight of the gravimetric water content at saturation $\rho_dgw_{sat}$ [kg.m⁻³.kg⁻¹.m³],

$$s = \frac{w}{w_{sat}}; \qquad (1.12)$$

* saturation $s$ is expressed as the ratio of the volumetric water content $\theta$ and the porosity $\phi$,

$$s = \frac{\theta}{\phi}. \qquad (1.13)$$

Note that the moisture index associated with saturation is not relevant for media whose porosity changes in the presence of water (e.g. clay media that swell in contact with water). In our study, we will assume that the medium porosity does not change in the presence of water, or that it is balanced on average in the medium.
These moisture indicators describe the distribution of water in a porous medium. There are three main classes of water present in a porous medium (Figure 1.2):

* gravitatory flow driven water is water flowing rapidly under gravitational forces, and is mainly found in macroporosity;

* the available water content is the water retained by the attractive force of the solids present mainly in the microporosity;

* capillary water is water strongly bound to solid particles.

<figure><img src="image_9.png" /><figcaption>Figure 1.2: Sketch of water distribution in a porous medium.</figcaption></figure>

In the case where the porous medium is a soil, this delimitation of water descriptions is associated with water accessibility by a plant. We call "field capacity" a soil without gravity water. Thus, we define the available water content as the water between the field capacity and the point of permanent wilting (plant death). The water in the available water content is available to the plant through root suction forces. Note that temporary wilting points correspond to transpiration by the plant exceeding water absorption in the soil (reversible wilting). In the case where water is not accessible by the roots, the water flows is rapidly into the macroporosity (gravity flow driven water) or it is retained by the solid phase (capillary water). This water is therefore unavailable to the plant.

A porous medium is never totally dry or totally saturated (all pores filled with water). A residual saturation remains adsorbed to the solid elements of the medium, and the medium is rarely saturated at 100%, as a residual amount of air is almost always trapped even in very wet pores [170]. We assume that the saturated water content $\theta_{sat}$ [m$^3$.m$^{-3}$] is approximately equal to the porosity, i.e. $\theta_{sat} \approx \phi$.

An effective moisture is then used to represent the water in the medium, excluding capillary water and the fraction of potentially trapped air flowing through the medium. This new indicator is generally expressed as a volume ratio. The effective volumetric water content $\theta_e$ [m$^3$.m$^{-3}$] is then expressed as the difference between the water content in the medium $\theta$ and the residual water content $\theta_r$,

$$ \theta_e = \theta - \theta_r. \qquad (1.14) $$

An effective saturation $s_e$ [m$^3$.m$^{-3}$] is generally expressed as the ratio between the effective volumetric water content and the difference between the water content at saturation $\theta_{sat}$ and the residual volumetric water content,

$$ s_e = \frac{\theta_e}{\theta_{sat} - \theta_r}. \qquad (1.15) $$
However, this indicator does not take into account the correct void volume, since it considers capillary water to be associated with the solid phase. For the remainder of this work, it is important to note that in multi-phase media such as porous media: the saturation of fluids in the pore space must be decomposed such that the sum of these saturation is 1. In our case, we consider water in the two forms described above water (bound or free) and air. Thus, the saturation of water s and the saturation of air s_a [-] are the only saturation in the pore space (Figure 1.3), i.e.

$$s + s_a = 1. \qquad (1.16)$$

The water saturation s in the medium is expressed as the sum of the residual water saturation

<figure><img src="image_10.png" /><figcaption>Figure 1.3: Sketch of a porous medium with different phase of fluid.</figcaption></figure>

$s_r$ [m$^3$.m$^{-3}$] and the free water saturation $s_w$ [m$^3$.m$^{-3}$],

$$s = s_r + s_w. \qquad (1.17)$$

In the presence of air in the pore space, water saturation has a limit which takes into account the air trapped in the medium. This limit is given by a compression threshold $s_{comp}$ such that $s_{comp} = 1 - s_a$. Note that if all the pores are filled with water, then the water content by volume at saturation is equivalent to the porosity, i.e. $\theta_{sat} = \phi$. Thus, water saturation admits a limit under air compression ($s_{comp} = 1$).

In the rest of this manuscript, we will mainly use the moisture indicator associated to the volume to facilitate modeling and define definition domains for our moisture in the porous medium.

### 1.1.2 Organic matter: a substance influencing soil structure

Soil is a porous medium composed of organic matter produced by living soil organisms. Figure 1.4 is a simplified sketch of a soil including the structure of a porous medium in the presence of organic matter as describe in Figures 1.5-1.6. The presence of organic matter and living organisms add complexity to the study of soil structure. In order to understand how organic matter interacts with the porous medium, it is essential to understand the living organisms and solid structures in the soil.

#### The living soil organisms

Soil is a porous medium colonized by a large variety of living organisms. They represent the entire pedoflora and pedofauna.
<figure><img src="image_11.png" /><figcaption>Figure 1.4: Sketch of a soil structure</figcaption></figure>

Soil fauna represents all the animals living in the soil. Depending on their size, soil fauna includes microfauna (organisms smaller than 0.2 mm, e.g. roundworms and protozoa), mesofauna (organisms between 0.2 and 4 mm, e.g. mites, springtails, enchytreids, nematodes), macrofauna (organisms between 4 and 80 mm, e.g. worms, insects, molluscs) and megafauna (organisms larger than 80 mm, e.g. moles, groundhogs, rabbits) [217, 55].

Pedoflora is the group of microorganisms and macroscopic underground plant organs living in the soil. Microorganisms (bacteria, fungi) make up the majority of living soil organisms [150]. They account for 1 to 4% of soil organic carbon [312]. Among them, fungi account for up to 80% of microbial biomass [132].

## Aggregates in the soil: the habitat of living organisms

Soil is a habitat produced by living soil organisms. It is made up of several heterogeneous assemblages of substances or elements that adhere strongly to each other, called aggregates. These structures are more or less large in the medium. They have complex and irregular organizational structure. They are classified by aggregation size [116, 330, 151]: large aggregates and lumps (greater than 2000 mm), macroaggregates (between 250 µm and 2000 mm), microaggregates (between 2 and 250 µm) and primary particles (less than 2 µm).

These aggregates offer liveable structures for soil organisms. They live in different soil areas in the void spaces adapted to their size. Therefore, the size of the living organisms present in a structure is correlated with the size of the aggregates. For example, microorganisms live in all structures larger than 0.2 µm, while megafauna and macrofauna live in large aggregates and lumps. Note that structures smaller than 0.2 µm are not habitable [217].

Pedofauna and pedoflora are important factors in soil aggregation. They constitute an important fraction of the aggregates and structure of their habitat. Mesofauna, macrofauna, megafauna and roots generate porosity to live in soil structures. For example, different species of earthworms burrow tunnels [184, 60] and roots move aggregates [278, 174, 153].

## Organic matter: a soil-aggregating substance

Living soil organisms make up and contribute to the living soil biomass. They produce organic matter, which accounts for 1 to 3% of soil solids [251]. Organic matter feeds a biotransformation cycle of living and dead soil biomass. It is in perpetual evolution, passing through three successive stages:
fresh organic matter (dead plants and animals), transitional or young humic products (mucilage, an Extracellular Polymeric Substances) and stable humus (cellulose and lignin). The biomass of living soil organisms represents around 15% of soil carbon [55]. Despite the low representation of living biomass, it plays an important role in the description of soil at different scales.

Organic matter produced by living soil organisms has a direct or indirect influence on soil structure. The release of viscous organic matter promotes the direct adhesion of substances to solid minerals and living or dead organisms forming aggregates. Microorganisms are the main agents aggregating the soil [78]. Among them, bacteria and organic substances such as mucilage form microaggregates [74] while fungi produce hyphae that hold aggregates together forming macroaggregates [94].

Roots also play an aggregating role. They liberate exudates, border cells [25, 133] and root mucilage [185] that affect aggregation formation and stabilize it around the root [278, 174, 153]. Root mucilage aggregates and lubricates the medium to promote root penetration into the structure. Root tips displace soil particles [220, 268], reducing mechanical impedance [185].

Other plant organic substances (exudates and root debris) are integrated with the microbial biomass to form biofilms. These biofilms contain Extracellular Polymeric Substances produced by microbial which act as binding agents in the microaggregates [76]. A similar effect is seen with pedofauna (earth-worms), which produce faecal matter and mucus that stimulate aggregation by the microbial biomass. This aggregation changes the microbial habitat into a more porous, ordered and aggregated structure [129]. The microaggregate stabilized by microbial matter, particularly polysaccharides, remains stable after the death of its hosts [333]. These processes are amplified by fine roots, which cause intense, periodic drying of their environment and promote aggregation processes [329, 278]. The drying phase is essential for aggregation. In the long term, this affects soil structure heterogeneity development in the soil profile [364]. In addition, interactions can accumulate during aggregate formation. For example, roots and associated mycorrhizal fungi promote soil macroaggregation [330, 94, 286, 353].

This overview of soil structure and the organic substances that affect it, proves that organic matter is a key substance in soil aggregation. In order to understand water transfer in the pore space, it is essential to look at the physical properties of this complex medium. A particular focus is set on the rhizosphere.

### 1.1.3 The hydraulic properties of a soil: application to the rhizosphere

The rhizosphere is the soil zone under the influence of plant roots [171, 174]. This zone is one of the most complex and crucial interfaces on Earth [174]. A multitude of biological, physical and chemical processes support all the other trophic levels of the biosphere [298].

#### Processes in the rhizosphere

Processes in the rhizosphere occur on multiple spatial and temporal scales. They result from biological, physical and chemical interactions between organisms, solids and fluids in the rhizosphere. These processes include, in particular: the development of the root hydraulic architecture [368], the impact of salinity stress [300, 187], the uptake of nutrients [114, 210], the interactions between soil structure and root growth [205, 220, 268], the impact of mucilage in root water uptake [301], the impact of organic matter on the soil structure [24], the behavior of organic matter in presence of microorganisms [254], the
uptake of pesticide [186] or the uptake of nutrients [114, 210].

Understanding processes in the rhizosphere should lead to better comprehend those at the woodland edge, forest edge and in the soil in general. For example, it is important to better understand the effects of organic matter on the structure and hydraulic properties of the medium.

## Organic matter in the rhizosphere

The rhizosphere is mainly composed of rhizodeposits, including root mucilage, which stimulates bacterial activity (degradation of organic matter and production of Extracellular Polymeric Substances) [265]. This mucilage is secreted at the root tip [185] and varies along root axis [179, 344]. The chemical composition of root mucilage can influence the hydraulic properties of the medium. This composition varies according to plant species [376, 245, 244, 37], microbial community [370] but also during their growth [72, 248, 58]. Carminati [62] has shown experimentally that mucilage-type substances have an impact on water content in the rhizosphere of 2 to 3 week-old roots. However, according to Nguyen *et al.* [249] and Mary *et al.* [229], the decomposition half-life of mucilage can vary between 3 and 11 days. This could be explained by the fact that bacteria degrade mucilage and produce substances similar to those of root mucilage [254].

## Hydraulic properties of mucilage in soil

Root mucilage is a gelatinous substance comprising mostly polysaccharides (between 78.4% [246] and 94% [38]). These polysaccharides are sugar units with amphiphilic properties, i.e. they have hydrophilic head (capillary effect) and hydrophobic tail (repellent water) which affect the hydraulic properties of the mucilage. This root mucilage impacts soil properties similarly as the Extracellular Polymeric Substances of bacteria [79, 367]. Here, we assume that we are not differentiating between root mucilage and the Extracellular Polymeric Substances of bacteria. These properties vary according to mucilage concentration and water saturation [277], i.e. more or less diluted mucilage changes its physico-chemical properties. These properties go from liquid to viscoelastic hydrogel [40], which influences the behavior of water in the soil and its structure [63, 361]. At high concentrations, when the medium is drying,

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

Figure 1.5: Examples of polymeric structures formed by mucilage and extracellular polymeric substances in porous media: (a) light microscope image of threads of mucilage (mucilage content = 4.5 mg.g⁻¹ [mg dry mucilage per g of particles]) formed across a large pore during drying; (b) light microscope image of cylinder formed between neighboring glass beads (1.7-2 mm in diameter) at intermediate mucilage content (0.7 mg.g⁻¹). Figure extract from Bernard *et al* [23].

the mucilage becomes solid and forms liquid bridges between the soil particles [23, 124, 125]. Figure 1.5
<figure><img src="image_13.png" /><figcaption>Figure 1.6: (a) Cutout of a horizontal slice of CT scan (chia pure) showing the filamentous network within air-filled voids including agglutinations (red arrows) of hydrated material. Larger hydraulic bridge indicated by blue arrow. Filamentous network was visualized by projecting the maximum gray value over 208 µm along the z-axis. Hydraulic bridge between chia seeds indicated by blue arrow. Pixel size in x- and y-direction is 20.8 µm. (b–d) Light microscope images of dehydrated chia mucilage in Nafion (b) and in sand (c,d). Figure extract from Van Veelen *et al* [344].</figcaption></figure>

and Figure 1.6 are schemes from Bernard *et al* [23] and Van Veelen *et al* [344] which represents liquid bridges due to mucilage in the porous medium, at pore scale. It illustrates the principle of aggregation to form a structure around the root, called rhizosheath. This aggregation is generated by the low surface tension and high viscosity of the mucilage [63]. This process is referred to as adsorption of the mucilage to soil solids. The mucilage is too sticky to diffuse into the pore space [277] due to the increased viscosity of the fluid [282, 130, 131, 314, 245] as saturation decreases. The glue effect of the mucilage forms a three-dimensional gel network in the pore space between the water and the solid [287, 131, 199, 37]. The forces in the soil increase [245, 252], improving connectivity between soil particles and water [66, 3, 63]. Because mucilage increases viscosity in soil, hydraulic conductivity decreases in the presence of mucilage for a constant load [238, 68, 199].

By increasing surface tension and decreasing viscosity, i.e the mucilage concentration decrease [281, 274], liquid bridges are broken [65, 68, 5, 298]. In other words, at low concentrations and in the presence of water, mucilage forms a Newtonian fluid with low viscosity [63, 124, 125, 298], which reduces the attraction with the solid of the liquid bridge. In this state, the mucilage promotes diffusion [277] that lubricates the soil and reduces friction between soil particles and the root, i.e. it facilitates root penetration into the soil.

Mucilage affects soil water flux during drying cycles that modify the water adsorption, viscosity, wettability [348]. It has a high water retention capacity [288, 280, 130, 302], swelling on contact with water and can hold water up to 1000 times its dry weigh [231]. It improves water storage compared with bulk soil during drying of the medium [201, 125]. However, it becomes hydrophobic during dry phase [61, 238, 24], i.e. repellent water occurs during dry cycle. Thus, the mucilage acts as a buffer, regulating hydraulic flow to the root. For example, rehumectation of the dry rhizosphere can take between a few hours and a few days [282, 65, 38, 244], illustrating the slow movement of water in this
rather hydrophobic medium.

This phenomenon is due to the fact that the hydrophobic tails of ExoPolySaccharides in an unsaturated medium are turned towards the gas phase. During humectation with this configuration, it is difficult for water to infiltrate through this hydrophobic zone. When water infiltrates, the exopolysaccharides configuration gradually changes, reducing the repulsive forces at the interface between the zone with and without exopolysaccharides (Figure 1.7).

In short, mucilage affects water transfer in the rhizosphere. To understand their influence on the rhizosphere processes, mathematical and numerical models of water transfer in their presence are essential.

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

Figure 1.7: Configuration of the liquid phase in soils containing Extracellular Polymeric Substances or mucilage. (a) In this illustration, the concentration of Extracellular Polymeric Substances or mucilage increases from the right to the left. During drying, the gas-liquid interface retreats and the polymers accumulate at this interface. At low polymer contents, the gas-liquid interface retreats but the liquid phase is not broken, which results in the formation of thin threads. At higher polymer contents, the gas-liquid interface becomes stiffer because of the entanglement of the polymers among themselves and with the soil particles. As an additional matrix. Together with the hygroscopic nature of the polymers, this leads to an amplified soil water retention. Besides increasing the water retention, the polymer network preserves the continuity of the liquid phase (the flow of water is illustrated by the dashed red arrows). (b) The liquid phase remains connected during drying, with the liquid converging into the two-dimensional surfaces. This induces a shift toward higher hydraulic conductivity in dry soils. Figure extract from Benard *et al* [23].
## 1.2 Model of water transfer

### 1.2.1 Model without amphiphilic matter

#### 1.2.1.1 Saturated medium: Darcy law

Darcy law [90] is a fundamental empirical relationship used to study a laminar, viscosity-dominated flow (low Reynolds number $Re \ll 1$) in a porous medium. It represents the water flow whose volume flow rate $Q$ [m³.s⁻¹] passes through a water-saturated porous medium of cross-section $S$ [m²] and length $L$ [m]. It establishes a proportional relationship between the volumetric flow rate $Q$ and the applied pressure gradient $P$,

$$Q = \frac{k}{\mu} S \frac{P}{L}, \qquad (1.18)$$

where $\mu$ is the dynamic viscosity of the flowing fluid [Pa.s] and $k$ is the intrinsic permeability of the medium [m²]. Figure 1.8 displays the Darcy experiment associated to the Darcy law for a constant potential $H$ (difference in water level between two reservoirs). Permeability is a variable parameter

<figure><img src="image_15.png" /><figcaption>Figure 1.8: Sketch of the Darcy experiment. Two columns containing water with a height differential $H$ are connected by a soil column of length $L$ and cross-section $S$. The system is assumed to include an open area for releasing a volume flow of water $Q$. This experiment measures the volume flow rate of the medium at a given pressure.</figcaption></figure>

that depends on the structure (grain size, grain shape and porosity). The hydraulic conductivity at saturation $K_{sat}$ defines the flow capacity of a fluid in a porous medium. This coefficient of permeability depends on the medium and the fluid passing through it:

$$K_{sat} = \frac{k \rho g}{\mu}, \qquad (1.19)$$

where $\rho$ is the density of the liquid [kg.m⁻³] and $g$ is the gravity [m.s⁻²]. Remember that hydraulic pressure $P$ is proportional to hydraulic head $H$,

$$P = \rho g H. \qquad (1.20)$$
Inserting Eqs (1.19)-(1.20) in Eq (1.18), Darcy law becomes

$$Q = K_{sat} \frac{SH}{L}. \qquad (1.21)$$

Darcy law is generalized to a velocity [167]:

$$v = -K_{sat} \nabla H, \qquad (1.22)$$

where $v$ is the mean velocity of the fluid, called Darcy velocity [m.s⁻¹].

Note that the actual velocity $v_r$ in the medium is the ratio between the mean velocity and the porosity of the medium:

$$v_r = \frac{v}{\phi}. \qquad (1.23)$$

This law is relevant in a homogeneous saturated medium for laminar flow with an incompressible Newtonian fluid. In other words, it is relevant as long as viscous forces dominate inertial forces.

### 1.2.1.2 Unsaturated porous medium

Buckingham [48] proposed a generalization of Darcy law for unsaturated media in 1907. He proposed that the hydraulic conductivity $K$ becomes a function of the water content of the porous medium. Thus, he established Buckingham law, also known as Buckingham/Darcy law, expressed:

$$v = -K(s) \nabla H. \qquad (1.24)$$

#### Hydraulic potential of the soil

Hydraulic potential is the energy state of water in the soil, which govern the movement of the liquid phase in the soil. This potential is affected by any potential difference that causes water to move in the soil. The total hydraulic potential of a soil, denoted $\Phi_t$ [J.kg⁻¹], corresponds to the sum of the components influencing the forces acting on the water and affecting its energetic state [260, 170, 256]. These components are generally:

* gravitational potential, denoted $\Phi_g$, corresponds to the gravitational forces.

* matric potential or capillary potential, denoted $\Phi_m$, corresponds to adsorption and capillary effects of the solid phase which interact with water phase. It is composed of capillary forces and attractive forces at the surface of solid particles. These forces are responsible for the retention of water in soil pores.

* air pressure potential, denoted $\Phi_a$, resulting from external gas pressure applied to the water. This potential is generally neglected.

* osmotic potential, denoted $\Phi_{\Pi}$, due to solutes in water. Molecules of water will move from the area with the highest concentration of solutes to the area with the lowest concentration, reducing the hydraulic potential. In soil, osmotic potential is often neglected except near root cell walls.
Therefore the total hydraulic potential is given by the sum of these components:

$$ \Phi_t = \Phi_g + \Phi_m + \Phi_{II} + \Phi_a + \dots \qquad (1.25) $$

Note that hydraulic potential is an energy that can be expressed in different ways [55]:

* hydraulic potential per unit mass $\Phi$ [J.kg⁻¹];

* hydraulic potential per unit volume of water $\Psi = \rho_w \Phi = \rho_w g \varphi$ [J.m⁻³]. This is homogeneous to a pressure [Pa];

* hydraulic potential per unit weight of water $\varphi = \frac{\Psi}{\rho_w g} = \frac{\Phi}{g}$ [J.N⁻¹]. This is homogeneous to a pressure head [m].

This concept of hydraulic potential in a soil is generally well defined in theoretical articles mentioning free energy [260, 170, 256]. However, there remains an ambiguity with measured pressure, which has the same physical dimension as hydraulic potential. A confusion arises where pressure is associated with hydraulic potential. We will discuss this in more detail in chapter 2 in order to clarify it and its components.

## Hydrodynamic properties of a porous medium

Water transfer in a porous medium is characterized by hydrodynamic properties. These properties are defined by the water retention characteristic curve $h(\theta)$ and the hydraulic conductivity $K(\theta)$ or $K(h)$.

**The water retention characteristic curve** is the representation of the soil's capacity to retain water. This retention is due to the water potential, also called matric potential, noted $h$ or $\psi_m$. This term is called potential here, but it corresponds to the pressure associated with the volumetric energy required to retain the water in the medium. We will discuss about this ambiguity further in the chapter 2. In this chapter, we assume that matric potential is a pressure.

In unsaturated medium, this matric potential is negative. Note that, the negative of this potential is called suction. As the soil becomes drier, the medium loses saturation, and the remaining water is more strongly bound to the soil. It becomes harder to extract it from the soil. More energy is required to loosen this water.

The retention curve is a measure of the bonding intensity between the water and solid phase. It varies for each soil with the characteristics of pore (size, shape and distribution) and decreases with the rate of water saturation. Then, the matric potential is related to water content by empirical function defined by the retention curve. In this study, we assume that it is expressed per unit weight of water, as is generally the case in the soil physics literature. In other words, it is expressed in terms of water column height.

Several empirical relationships have been established over the last few decades, as shown in Leong's, Bullied's or Too's reviews [211, 49, 332]. These functions have to be fitted with retention measurements, i.e some parameters need to be fixed. The number of parameters varies according to the model chosen (from 2 parameters for the Williams model to 5 parameters for the Omuto model). The table 1.1 summarizes the various relationships associated with the water retention characteristic curve and the
number of parameters to be fixed. Among them, the Brooks and Corey model and van Genuchten model are the most used. Note that the parameter *m* in the van Genuchten model is related to the parameter *n*. Mualem [240] and Burdine [50] propose to relate these two parameters through relationships between the retention curve and hydraulic conductivity. As a result, relationships between these two parameters have been established:

$$m = \begin{cases} 1 - \frac{1}{n}, & \text{Mualem,} \\ 1 - \frac{2}{n}, & \text{Burdine.} \end{cases} \qquad (1.26)$$

<table><caption>Table 1.1: Common water retention characteristic models in literature</caption><thead><tr><th>Model</th><th>Function θ(h)</th><th>Parameters</th><th>Remarks</th></tr></thead><tbody><tr><td>Gardner (1958)</td><td>θ<sub>r</sub> + (θ<sub>sat</sub> - θ<sub>r</sub>) (1 + αh<sup>n</sup>)<sup>-1</sup></td><td>θ<sub>sat</sub>, θ<sub>r</sub>, n, α</td><td>[145]</td></tr><tr><td>Brooks and Corey (1964)</td><td>θ<sub>r</sub> + (θ<sub>sat</sub> - θ<sub>r</sub>) (H<sub>cap</sub>/h)<sup>λ</sup></td><td>θ<sub>sat</sub>, θ<sub>r</sub>, λ, H<sub>cap</sub></td><td>[43]</td></tr><tr><td>Farrel and Larson (1972)</td><td>1 - (1/α) ln(h/a)</td><td>α, a</td><td>[126]</td></tr><tr><td>Campbell (1974)</td><td>θ<sub>sat</sub> (h<sub>b</sub>/h)<sup>1/b</sup></td><td>θ<sub>sat</sub>, h<sub>b</sub>, b</td><td>[57]</td></tr><tr><td>Van Genuchten (1980)</td><td>θ<sub>r</sub> + (θ<sub>sat</sub> - θ<sub>r</sub>) [1 + (αh)<sup>n</sup>]<sup>-m</sup></td><td>θ<sub>sat</sub>, θ<sub>r</sub>, n, m, α</td><td>[343, 240, 50]</td></tr><tr><td>Tani (1982)</td><td>θ<sub>r</sub> + (θ<sub>sat</sub> - θ<sub>r</sub>) [1 + αh exp(-αh)]</td><td>θ<sub>sat</sub>, θ<sub>r</sub>, α</td><td>[319]</td></tr><tr><td>Williams (1983)</td><td>exp[ln(h)/b]</td><td>a, b</td><td>[360]</td></tr><tr><td>McKee and Bumb (1984)</td><td>θ<sub>r</sub> + (θ<sub>sat</sub> - θ<sub>r</sub>) exp(a/h)</td><td>θ<sub>sat</sub>, θ<sub>r</sub>, a, b</td><td>[233]</td></tr><tr><td>McKee and Bumb (1987)</td><td>θ<sub>r</sub> + (θ<sub>sat</sub> - θ<sub>r</sub>) [1 + exp(a/h)]<sup>-1</sup></td><td>θ<sub>sat</sub>, θ<sub>r</sub>, a, b</td><td>[232]</td></tr><tr><td>Russo (1988)</td><td>θ<sub>r</sub> + (θ<sub>sat</sub> - θ<sub>r</sub>) [(1 + αh/2) exp(a/h/2)]<sup>2/(m+2)</sup></td><td>θ<sub>sat</sub>, θ<sub>r</sub>, α, m</td><td>[295]</td></tr><tr><td>Fredlung and Xing (1994)</td><td>θ<sub>r</sub> + (θ<sub>sat</sub> - θ<sub>r</sub>) [ln(e + (h/a)<sup>n</sup>)]<sup>-m</sup></td><td>θ<sub>sat</sub>, θ<sub>r</sub>, a, m</td><td>[136]</td></tr><tr><td>Kosugi (1999)</td><td>θ<sub>r</sub> + (θ<sub>sat</sub> - θ<sub>r</sub>) erfc[ln(h)/σ√2]</td><td>θ<sub>sat</sub>, θ<sub>r</sub>, h<sub>m</sub>, σ</td><td>[196]</td></tr><tr><td>Groenevelt and Grant (2004)</td><td>θ<sub>r</sub> + k<sub>1</sub> [exp(-k<sub>0</sub>/15000<sup>n</sup>) - exp(-k<sub>0</sub>/h<sup>n</sup>)]</td><td>θ<sub>r</sub>, k<sub>0</sub>, k<sub>1</sub>, n</td><td>[156]</td></tr><tr><td>Omuto (2009)</td><td>θ<sub>r</sub> + θ<sub>1</sub> exp(-α<sub>1</sub>h) + θ<sub>2</sub> exp(-α<sub>2</sub>h)</td><td>θ<sub>r</sub>, θ<sub>1</sub>, θ<sub>2</sub>, α<sub>1</sub>, α<sub>2</sub></td><td>[253]</td></tr></tbody></table>

Despite the large number of models, all the h(θ) curves are decreasing, so for a given h matric potential, there is a single h(θ) water content. However, experiments show an hysteresis phenomenon during the wetting/drying cycles of the medium. Figure 1.9 illustrates schematically the hysteresis of this retention curve. There are different paths representing retention as a function of water content. The curves in the drying phase are shown in red, while the curves associated to the wetting phase are shown in blue. The solid lines represent the main curves of each phase, while the dotted lines represent the transient curves during a changing phase during a cycle. This phenomenon, which occurs over
<figure><img src="image_16.png" /><figcaption>Figure 1.9: Sketch of a water retention curve during a wetting/drying cycle.</figcaption></figure>

several cycles, makes the study complex.

The impact of organic matter on these retention curves and the hysteresis phenomenon has not yet been fully modelled.

**The hydraulic conductivity** is the capacity of a porous medium to infiltrate water, and relates the water content to the hydraulic conductivity $K$. This expresses a soil's ability to let a fluid through under the effect of a pressure gradient.

The hydraulic conductivity is the ratio between the volume flow of water (expressed in m.s⁻¹) through a given soil section and the gradient head (expressed without unit) governing this flow: Darcy law if the medium is saturated, or Buckingham's law if the medium is unsaturated. Therefore, it is homogeneous to a water flow or velocity, and is expressed in m.s⁻¹.

In a saturated medium, all pores are filled and hydraulic conductivity is at its maximum ($K = K_{sat}$). By analogy, the pore network can be assimilated to a set of capillaries, which allow water to flow through them. However, the forces of the substrate become more and more difficult to transfer water when the soil is drying. So, during drying processes, the pores with the largest diameters are emptied first: the water flows into pores with smaller and smaller diameters. As a result, friction forces increase, slowing down the flow through smaller pores and reducing hydraulic conductivity. Water flow is much lower in dry soil than in moist soil. Hydraulic conductivity is a function proportional to hydraulic conductivity, and relative conductivity $k_r$ (increases with soil water saturation),

$$K(s) = K_{sat}k_r(s). \qquad (1.27)$$

This model most commonly used to describe the hydraulic conductivity curve is that proposed by Van Genuchten [343] using Mualem model [240],

$$K(s) = K_{sat}s^l \left(1 - \left[1 - s^{\frac{1}{m}}\right]^m\right)^2, \qquad (1.28)$$

where $K_{sat}$ is the hydraulic conductivity at saturation [m.s⁻¹] and $l, m$ are parameters. Note that generally $l = 0.5$ and $m = 1 - \frac{1}{n}$ [240]. Figure 1.10 illustrates the hydraulic conductivity curve for a
saturation hydraulic conductivity of $10^{-4}$ for different shape parameters $n$.

<figure><img src="image_17.png" /><figcaption>Figure 1.10: Example of hydraulic conductivity using Van Genuchten model (Eq. (1.28)). Parameters $K_{sat} = 10^{-4}$ m.s$^{-1}$, $l = 0.5$ and $m = 1 - \frac{1}{n}$.</figcaption></figure>

Other relationships are available in the literature, such as the Brooks and Corey model [44] or the Burdine model [50]. Nevertheless, all the conductivity curves $K(s)$ are an increasing function and do not differ qualitatively.

Richards equation

Water transfer in an unsaturated porous medium is described by the Richards equation. In 1931, Richards [284] proposed coupling the Buckingham-Darcy equation with the mass conservation equation to describe dynamic phenomena in a porous medium. The mass conservation equation for a liquid is given by

$$ \partial_t \rho + \nabla \cdot \rho \vec{U} = 0, \qquad (1.29) $$

with $\rho$ the density of the fluid and $\mathbf{U}$ the velocity of a fluid particle. The density of the fluid in a porous medium $\rho_p$ depends on the density of the water $\rho_w$, the porosity $\phi$ and the saturation $s$:

$$ \rho_p = \rho_w s \phi. \qquad (1.30) $$

Thus, the mass conservation of the water phase in the porous medium reads

$$ \partial_t (\phi \rho_w s) + \nabla \cdot (\rho_w \mathbf{v}) = 0. \qquad (1.31) $$

The medium is assumed to not be deformed, or to be slowly deformed compared to the fluid motions, i.e. the porosity is assumed to be constant. Inserting Darcy law Eq. (1.18) into the mass conservation of the water phase Eq. (1.31), we obtain Richards equation

$$ \left\{ \begin{array}{l} \phi \partial_t s + \nabla \cdot [K(s) \nabla H] = 0, \\ H = h(s) + z. \end{array} \right. \qquad (1.32) $$

Note that this formulation is equivalent to the following:

$$ \partial_t \theta + \nabla \cdot [K(\theta) \nabla (h_c(\theta) + z)] = 0. \qquad (1.33) $$
In soil physics, the hydraulic head associated with capillary pressure is defined by the water retention curve. This curve is monotonic and decreasing. A bijective function exists between water content and hydraulic head $h_c$. Using this relationship between matric potential and water content, Richards equation can be written as follows

$$C(\theta)\partial_t h_c + \nabla \cdot [K(h_c)\nabla(h_c + z)] = 0, \qquad (1.34)$$

with $C(\theta)$ specific capillary capacity [m⁻¹]

$$C(\theta) = \frac{\partial\theta}{\partial h}. \qquad (1.35)$$

Note that Richards equation can also be expressed as

$$\partial_t \theta + \nabla \cdot [D(\theta)\nabla \theta] + \nabla \cdot [K(\theta)\vec{e}_z] = 0, \qquad (1.36)$$

with $D(\theta)$ the specific diffusion coefficient [m².s⁻¹]

$$D(\theta) = K(\theta) \frac{\partial h}{\partial \theta}. \qquad (1.37)$$

This equation is a complex relationship in the natural sciences [127]. Despite its similarity to the diffusion equation, its non-linearity and variations in diffusion according to water content make it complex to study. There are methods for solving this equation on one-dimensional domains and obtaining analytic solutions [345]: Kirchhoff transform, Laplace transform, Boltzmann transform or Traveling wave solutions. These methods have been also used to solve higher-dimensional problems using one-dimensional formulations. This model is one of the most widely used to represent transfer in porous media, despite its poor representation in highly unsaturated media. Note that there are other water transfer models based on other laws of fluid mechanics in porous media, such as the Navier-Stokes, Stokes or Brinkman equations (high porosity $\phi > 0.95$ [115]).

We have seen many models have been developed to represent water transfer in porous media. Therefore, we can investigate the development of water transfer models that take into account the presence of organic matter (mucilage), which influences the transfer (1.1.3).

### 1.2.2 Model in effective medium in presence of mucilage

One of the first approaches to modeling water transfer in the presence of mucilage is to consider a homogeneous medium and measure its water properties. For example, models have been developed for soil sampling that neglect the influence of mucilage [290, 183, 111]. Thus, they consider that transfer in this soil zone is modeled by the Richards model (section 1.2.1.2).

Mucilage introduces heterogeneity into the medium, particularly at the interfaces between soil with and without mucilage. These interfaces are generally considered to be at non-equilibrium due to swelling and shrinkage of the mucilage during wetting and drying cycles. Carminati [61] proposes adding the first order time constant equilibration model with Richards model to represent the wetting rate of mucilage in this non-equilibrium zone. He considers that the water content $\theta$ in the rhizosphere
is formed by a fraction of water content at equilibrium without mucilage $\theta_1$ and another fraction of water content with mucilage $\theta_2$ which is non-equilibrium i.e.

$$ \theta = \theta_1 + \theta_2 \qquad (1.38) $$

Note that at equilibrium, according to the retention curve, there is a unique matric potential $h_1$ associated with the water content $\theta_1$. Thus, he considers that the dynamics in the rhizosphere are described by

$$ \partial_t \theta = C(\theta_1)\partial_t h_1 + \Gamma(\theta) [h - h_{eq}], \qquad (1.39) $$

where $\theta$ is the current water content in the rhizosphere, $C$ the specific capillary capacity in the rhizosphere, $h$ the current matric potential in the rhizosphere, $h_{eq}$ the theoretical equilibrium matric potential for $\theta_2$ and $\Gamma$ is a proportional coefficient that describes the wetting rate of the rhizosphere.

This model was designed in analogy with the model used by Ross and Smettem [293] to model water content at the interface between the heterogeneities of aggregates and their sizes. The first-order time constant equilibration model was developed in particular for mobile and immobile solute transport in a soil zone [261] and between fracture macroporosity and matrix pores [147].

Carminati model takes into account a fraction of soil with mucilage (non-equilibrium dynamics) and one without, following the classical transfer equations, with a clear boundary between soil without and with mucilage. This model has been made more flexible and complete by Kroener et al. [201]. They consider a progressive transition between bulk soil and rhizosphere. They define a nonmucilage-filled pore ratio $R$ in the rhizosphere. This ratio is the ratio between the volume of pores that are not filled with mucilage and the total pore volume. Assuming that the rhizosphere is located between two extremes: the soil completely saturated with mucilage ($R = 0$) and without mucilage ($R = 1$). The water content in the rhizosphere is composed of the fraction of water content in soil without mucilage (bulk soil) $\theta_{bulk}$ and the fraction of water content with mucilage $\theta_M$,

$$ \theta = R\theta_{bulk} + (1 - R)\theta_M. \qquad (1.40) $$

They build a model in which the flow dynamics in the rhizosphere is an average between the dynamics with and without mucilage:

$$ \partial_t \theta = RC_{bulk}(h)\partial_t h + (1 - R) \frac{1}{\tau} (h - h_{eq}(\theta_M)) \qquad (1.41) $$

where $\theta_M$ is the volumetric water content in the fraction of the elementary volume ($1 - R$) that is influenced by mucilage. Note that $\Gamma = \frac{1}{\tau} = \tau_0\theta_M^{-\gamma}$ where $\tau_0$ and $\gamma$ are parameters fitted from empirical data.

Because mucilage has a water-holding capacity, it increases water storage. In other words, it influences the retention curve. Mucilage decreases the water potential for a given volume of water. This decrease is due to osmotic forces. They assume that the osmotic potential in the gel is related to the gravimetric mucilage concentration in water $c$ [g.g$^{-1}$] by a power law:

$$ \pi(c) = \omega_0 c^{\beta}, \qquad (1.42) $$
with the gravimetric concentration of mucilage in water is the ratio of the dry weight of mucilage to the weight of water in the mucilage-filled pore space,

$$c = \frac{c_{tot} \rho_{bulk}}{\theta_{M} \rho_{w}}. \qquad (1.43)$$

Thus, they model the effects of water retention for a concentration of mucilage in the soil with the following relationship:

$$h_M(\theta_M) = h(\theta) - w_0 c^\beta, \qquad (1.44)$$

where $h_M(\theta_M)$ and $h(\theta)$ are the pressure heads with and without the impact of mucilage, $\theta_M$ and $\theta$ are the volumetric water contents with and without the impact of mucilage and $w_0$ and $\beta$ are fitting parameters that describe how $c$ affects the overall water retention behavior.

The retention curve associated with gel is a modified version of the retention curve in soil without mucilage. This model was built to take into account the effects of mucilage during drying: increasing concentration leads to increased osmotic and capillary forces.

In parallel, Ahmed et al [3] have developed a model of hydraulic conductivity in the presence of mucilage. They admit that in the presence of polysaccharides in the liquid phase, the viscosity of the gel increases in the soil. In other words, the higher the mucilage concentration, the more viscous the gel is. As a result, hydraulic conductivity decreases in presence of polysaccharides:

$$K_M(\theta, \mu) = \frac{\mu_w}{\mu} K_{bulk}(\theta), \qquad (1.45)$$

where $K_M$ and $K_{bulk}$ are the hydraulic conductivities with and without mucilage, and $\mu$ and $\mu_w$ are the viscosities of the liquid phase with and without mucilage respectively. Note that the viscosity of the liquid phase without mucilage corresponds to water ($\mu_w = 1 \text{ mPa.s}^{-1}$). The viscosity of the liquid phase with mucilage is a control factor similar to the viscosity law of a polymer solution established by Morris et al [239] such that

$$\mu(c) = \mu_w (1 + v c^d), \qquad (1.46)$$

where $v$ and $d$ are two fitted parameters.

This hydraulic conductivity relationship was introduced into the Schwartz model [301] to improve the description of the Kroener model (Eq. (1.41)) for water uptake of roots. This model was made more flexible by Kroener et al [200], by introducing a fraction of soil with mucilage with a linearly decreasing concentration $R(c) = \exp(-ac)$, with $a$ a parameter to be adjusted. Mucilage concentration decreases from the root surface toward the bulk soil [214, 143].

Schwartz model is currently present in coupling models that take into account processes interacting with water transfer in root uptake, such as the R-SWMS model [183] or more recently in the model developed by Landl et al [206]. The latter model represents the water uptake of the root during elongation and introduces a mucilage concentration gradient in the rhizosphere.

Therefore, there are models that take into account the effects of mucilage in the rhizosphere for constant concentrations [61, 64, 301] and for variable concentrations between rhizosphere and bulk soil [200]. Because mucilage impacts soil parameters according to their spatial distribution, smaller-scale models have been developed to better model transfer in the presence of mucilage.
1.2.3 Pore scale model

During the last decade, several models have been developed, seeking to integrate biological, chemical
and physical processes in the rhizosphere [298]. Recent work has focused on better modeling the effects
of mucilage, particularly at the pore scale.

Jahromi *et al* [125] have developed a model based on a lattice spring method (LSM) to simulate the distribution of mucilage in the dry state at the pore scale. However, for wetter states, it is necessary to take into account additional physical phenomena such as surface tension, contact angle and viscoelasticity. They aim to develop a Lattice-Boltzmann simulation framework (methods commonly used to simulate pore scale dynamics of liquids [272, 285, 315, 340]) to simulate two-phase flow involving mucilage. This model should improve the representation of the retention curve in the presence of mucilage [125].

However, this type of model requires the geometry of the soil structure in order to be used. Imaging
methods (X-ray computed tomography [334, 275] or magnetic resonance imaging [313, 157]) are gener-
ally used to characterize the spatial distribution of soil structure (macropore networks and plant roots
in 3D) and to visualize and quantify soil water dynamics [297, 218]. Note that this question of imaging
coupled with predictive modeling is a recent topic for understanding rhizosphere processes [273, 291],
the challenges of modeling soil processes [346], or plant-soil modeling [294]. In general, current models
do not take into account the modification of geometry over time due to the actions of mucilage and
their degradation by microbial action.

Rötzer *et al* [296] have proposed a model representing soil aggregation dynamics in the presence of mucilage and a root. This model is an interesting perceptive for studying the evolution of soil structure, and could be coupled with one of the existing pore scale models. However, to reproduce the transfer dynamics on a larger scale, the challenge is to select a correct pore volume representative of the medium. A mathematical homogenization of the domain is required [279], but this is not always apparent. Thus defining this representative volume is easier in controlled media (microfluidic model [289, 311, 310]) but remains very complex to define in a natural medium such as the rhizosphere.

Pore scale models are very useful for better understanding certain local pore scale processes, but they are difficult to extend to larger scales. Moreover, models in the rhizosphere do not take into account the hydrophobic interactions present in the concept of soil wettability in the presence of amphiphilic or hydrophobic matter. Therefore, wettability offers an potential improvement to these water transfer models in the presence of amphiphilic matter.

1.3 Modeling of wettability and application to the soil

1.3.1 Wettability on plane substrate

Wetting is the study of the spreading of a liquid placed on a solid (or liquid) substrate [93]. Wettability is the ability of a fluid to spread on a surface. This phenomenon affects a wide range of fields, including industry, the automotive industry, pharmaceuticals and soil.

Understanding surface wetting means understanding why water spreads on hydrophilic substrates rather than hydrophobic surfaces. Moreover, it allows us to understand the spreading behavior of a
liquid in the presence of amphiphilic substances affecting the wettability of a surface. These effects are due to intermolecular forces resulting from capillary phenomena.

## Capillary phenomena

Capillary phenomena take many forms: the rise of water with nutrients by the roots in the stems of a plant, the formation of a meniscus in a glass of water or the rise of water in a sponge or wall [93]. They occur at interfaces between two non-miscible liquids, between a liquid and air, or between a liquid and a solid surface. They are due to surface tension forces corresponding to surface energy.

Surface tension is a physico-chemical phenomenon related to molecular interactions in a fluid. Molecules interact with each other in a fluid. However, in the neighbourhood of another phase (solid, liquid or gas), the molecules in this interface are influenced by the molecular forces of the other phase (Figure 1.11). The forces that interact are generally attractive (Van der Waals) and long-range. They can also be repulsive, but act at short range [227]. These interactions at this interface explain why a

<figure><img src="image_18.png" /><figcaption>Figure 1.11: Sketch of surface tension. The black arrows correspond to attractive molecular interactions.</figcaption></figure>

Gerris can stand on water and move along its surface without going underwater, or why a needle can float.

From an energetic point of view, surface tension influences the shape of the liquid. When a liquid is influenced by elementary work W [J] at an interface A [m²], it modifies its shape to adjust its surface tension and minimize its free surface (available surface not subject to viscous stress from the other phase). This elementary work is proportional to surface tension γ [N.m⁻¹],

$$dW = \gamma dA.$$

Thus, when the liquid supplies (resp. receives) elementary work, the molecules on the liquid surface have lost (resp. gained) their attractive interactions, affecting the surface tension of the liquid. This change represents an energy cost in order to maintain the liquid in its current state. Indeed, the surface tension at the γ interface corresponds to the derivative of the free energy of a closed system (liquid
interface in contact with an non-miscible phase) $\mathcal{F}$ with respect to the surface of the A interface:

$$ \gamma = \left[ \frac{\partial \mathcal{F}}{\partial A} \right]_{T,V,n}, $$

for a constant volume $V$ and quantity of matter $n$ at thermodynamic equilibrium $T$ [93]. Therefore, a liquid with a free surface minimizes its surface area until its surface tension minimizes the energy of the system. This explains why a drop with a free surface has a spherical shape. Thus, to increase the spherical surface area of the drop, in our example, we need to supply it with energy.

However, the curvature of the free surface implies the formation of over pressure and under-pressure on the liquid surface (pressure discontinuity), called Laplace pressure. This pressure tends to increase (in the case of over-pressure) or decrease (in the case of under-pressure) the water level at the convex or concave part respectively of the water surface (Figure 1.12) [19].

<figure><img src="image_19.png" /><figcaption>Figure 1.12: Sketch of Laplace pressure on the free surface (red line): over-pressure $p_+$ if the surface is convex and under-pressure $p_-$ if the surface is concave.</figcaption></figure>

Laplace pressure at the interface between a liquid and another non-miscible phase expresses the pressure difference between the inside of any fluid $p_{int}$ and the outside pressure $p_{out}$ at the interface as proportional to the mean surface curvature $\kappa$ and surface tension [93]:

$$ \Delta p = p_{int} - p_{out} = \kappa\gamma. $$

In the case of a liquid on a solid, the free surface of the liquid is curved on contact with a solid wall along a typical length, noted $l_c$: this length is called the capillary length. The curvature of the free liquid surface is governed by capillary forces and gravity. Thus, capillary length results from the equilibrium between Laplace pressure and hydrostatic pressure ($\Delta p = \rho g l_c$):

$$ l_c = \sqrt{\frac{\gamma}{\rho g}}. \qquad (1.47) $$

Note that the capillary length corresponds to the inverse of the curvature ($\kappa^{-1}$).

This length is the critical equilibrium length between gravity and capillary forces. When the fluid length is greater than the capillary length, gravity dominates capillary forces. Therefore, on a solid, flat, horizontal surface, a drop of radius $R \gg l_c$ is flattened by the gravity of its own weight. At
equilibrium, it adopts a puddle shape of thickness $e = 2\ell_c \sin(\frac{\Theta}{2})$ far from the triple point and forms a contact angle $\Theta$ with the substrate. Conversely, if its radius $R \ll \ell_c$, the 2D drop (reps. 3D) forms a parabolic cap (reps. spherical) which forms a contact angle $\Theta$ with the substrate. The notions of wetting presented in the next section, such as contact angle and triple line <fnref n="1" />, appear naturally here. These notions of wetting and capillarity are important for understanding the shape taken by the free surface of a liquid in contact with a solid. This is particularly true for the formation of a meniscus on the liquid surface between two plates of a vertical solid, or the rise of water in a tube with a small radius and a hydrophilic wall (Jurin's law).

Therefore, surface tension and capillary forces between a fluid and another non-miscible phase influence the shape of the liquid interface. In the case of solid substrates, this results in wetting phenomena such as spreading of the liquid on the solid substrate, or the formation of a contact angle.

## Spreading parameter

The spreading parameter $S$ characterizes the spreading of a liquid on a surface. It expresses the difference between the surface energies of the dry substrate $E_{\text{dry}}^{\text{substrate}}$ and the wet substrate $E_{\text{wet}}^{\text{substrate}}$,

$$S = E_{\text{dry}}^{\text{substrate}} - E_{\text{wet}}^{\text{substrate}}. \qquad (1.48)$$

If a drop is placed on a flat, smooth, homogeneous surface, it adapts its shape and spread to minimize its energy. This minimization results from competition between the interactions on each free surface: the surface tensions between the solid and the liquid $\gamma_{SL}$, the surface tensions between the solid and the gas $\gamma_{SG}$ and the surface tensions between the liquid and the gas $\gamma_{LG}$ or $\gamma$. Thus, the spreading parameter described in Eq. (1.48), can be rewritten,

$$S = \gamma_{SG} - (\gamma_{SL} + \gamma). \qquad (1.49)$$

This parameter expresses the spreading behavior of the liquid on the substrate and defines whether the substrate is wetting or non-wetting. When a drop is placed on a flat and hydrophilic surface, the drop spreads on the surface. Two types of spreading behaviour are distinguished:

* the drop spreads infinitely until it completely covers the surface, i.e. the surface energies of the dry substrate are dominant over the surface energies of the wet substrate, so the spreading parameter is positive ($S > 0$);

* the drop has a finite spread that does not cover the entire surface, i.e. the surface energies of the wetted substrate are dominant over the surface energies of the dry substrate, so the spread parameter is negative ($S < 0$).

Note that a liquid will spread completely if it is less polarizable than a solid. These two wetting regimes, represented by figure 1.13 are called complete wetting ($S > 0$) and partial wetting ($S < 0$).

However, the spreading parameter is not sufficient to highlight the behavior of the liquid on non-wetting surfaces. On hydrophobic surfaces, the liquid withdrawal until the surface area of the wetted substrate is almost vanished. In this case, the substrate is considered as non-wetting.

<footnote n="1">¹The triple line is an imaginary line resulting from the interactions of the three phases at the wettability front.</footnote>
<figure><img src="image_20.png" /><figcaption>Figure 1.13: Sketch of the complete and partial wetting depending to the spreading parameter S.</figcaption></figure>

Note that the theoretical criteria associated with the spreading parameter are related to the empirical criteria of Zisman [207, 135] who considers that any solid substrate has a critical surface tension $\gamma_c$: partial wetting ($\gamma > \gamma_c$) or complete ($\gamma < \gamma_c$). Therefore, considering that $\gamma_c \sim \gamma_{SG} - \gamma_{SL}$, the definition of wetting given by De Gennes [92] is found: complete wetting ($\gamma + \gamma_{SL} < \gamma_{SG}$) and partial wetting ($\gamma_{SL} < \gamma_{SG} < \gamma_{SL} + \gamma$). The medium is wettable if $\gamma_{SG} > \gamma_{SL}$ else it is not wettable ($\gamma_{SL} \ge \gamma_{SG}$).

However, this notion alone does not characterize all types of partial wetting, such as finite withdrawal (see Figure 1.14b). To represent these partial wettings, we need to define additional concepts.

### Triple line and contact angle

The triple line is an imaginary line resulting from the interactions of the three phases at the wettability front. At the junction of these surface interactions, a $\Theta$ angle is generally defined as a contact angle. This angle represents the degree of wettability of the liquid on the substrate. Note that this triple line and contact angle evolve during the drop spreads on the flat substrate.

At equilibrium, Young [365, 318] defined a relation between the wettability front and the static contact angle $\Theta_e$,

$$ \gamma \cos \Theta_e = \gamma_{SG} - \gamma_{SL}. \qquad (1.50) $$

Young relation can be rewritten as follows in order to relate the static contact angle and the spreading parameter,

$$ S = \gamma (\cos \Theta_e - 1). \qquad (1.51) $$

<figure><img src="image_21.png" /><figcaption>Figure 1.14: Sketch of partial wetting (a) hydrophilic and (b) hydrophobic tendency.</figcaption></figure>

This new concept is used to study the behavior of partial wetting. Figure 1.14a represents partial
wetting with hydrophilic (spreading) and hydrophobic (withdrawal) tendencies. A partial wetting is considered to have a hydrophilic tendency if the static contact angle satisfies the condition: $0 < \Theta_e < \frac{\pi}{2}$. Conversely, it is considered hydrophobic if $\frac{\pi}{2} \le \Theta_e < \pi$.

The notion of contact angle and spreading parameter are generally used for experimental characterization of substrate wettability, according to the following classification:

* hydrophilic substrate with complete wetting if $S > 0$ and $\Theta_e = 0$;

* hydrophilic substrate with partial wetting if $S < 0$ and $0 < \Theta_e < \frac{\pi}{2}$;

* hydrophobic substrate with partial wetting if $S < 0$ and $\frac{\pi}{2} \le \Theta_e < \pi$;

* non-wetting substrate if $S \le 0$ and $\Theta_e = \pi$.

### 1.3.2 Wettability in the soil

A porous medium has no continuous surface, unlike a flat substrate. Its surface is at least rough, and its geometry is highly complex. Wetting may therefore vary from the basic description given in section 1.3.1. De Gennes et al. [93] consider the "wetting of porous media" as being similar to the wetting of highly irregular surfaces, but with the component of imbibition. Thus a drop placed on a porous surface spreads and penetrates its support, whose modifies its wetting properties. A film deposited on a porous substrate is removed by two mechanisms: suction and dewetting. This domain is very large, covering the wetting of powders, rocks and soils, as well as the entire paper industry.

#### Capillary phenomena in soil

Porous media are formed by randomly distributed pores with complex geometries. They are often simplified by vertical capillaries of simple cross-section. In this case, the capillarity phenomena presented in section 1.3.1 occur. The phenomenon of capillarity is also present at the pore scale: wetting on a grain of sand, for example.

On macroscopic scales, these phenomena are less evident without this ideal framework. In petroleum engineering literature, the Leverett J-function was introduced by Leverett [212] to describe the capillary pressure $P_c$ [Pa] in a porous medium:

$$J(s) = \frac{P_c(s)}{\gamma \cos \Theta} \sqrt{\frac{k}{\phi}}. \qquad (1.52)$$

De Gennes et al. [93] proved, by macroscopic reasoning, that the capillary height $H_{cap}$ is defined in a porous medium such that

$$H_{cap} \approx \frac{\gamma \cos(\Theta_e)}{\rho g r_0} = \frac{\ell_c^2}{r_0} \cos(\Theta_e), \qquad (1.53)$$

where $\ell_c$ is the capillary length, $r_0$ is the typical pore radius such that $r_0 \approx \sqrt{\frac{k}{\phi}}$ and $\Theta_e$ is the equilibrium contact angle of the liquid on the non-porous substrate. Thus, capillary rise in a porous medium is similar to Jurin law [189], but expressed with macroscopic parameters of the porous medium [212, 93].
Note that the Leverett J-function and capillary height are linked by the relationship

$$P_c(s) = \rho g H_{cap} J(s). \qquad (1.54)$$

In addition, this notation is consistent with the models defined for represented water retention characteristic functions described in the Table 1.1 (section 1.2.1.2).

## Characterization of soil wettability

The wettability measurement of a porous medium is generally based on the direct measure of the contact angle at the surface of the porous medium or at the pore scale on an aggregate. This measure is used to characterize the capillary pressure.

At the pore and aggregate scale, measurement requires a tomography device and a representative pore volume. This type of measurement is useful in controlled microfluidic media [289, 311, 310] but is more difficult to define a representative pore volume in a natural environment.

On a macroscopic scale, the contact angle is measured at the surface of the porous medium. The wettability of the porous medium is defined at its surface. Note also that the effective contact angle at the surface [92] is not always the same as that at the pore scale [303]. In addition, the contact angle value is affected by hysteresis phenomenon, surface irregularities and the imbibition of the droplet in the medium. Particular wettability models are required to understand these phenomena: wetting models on physically textured surfaces (Wenzel model [358]) or chemically textured surfaces (Cassis-Baxter model [69]). These models are used to describe drop infiltration at the pore scale [35].

In the literature on powdery media, an alternative method is used to remove the imbibition component. This involves applying a thin layer of substrate to the tape and measuring the contact angle. However, this method remains contestable: the wettability of the tape can influence the contact angle measured.

Indirect measurement methods were developed to measure the contact angle at the surface of the medium: Molarity of Ethanol Droplet (MED) test [357, 194] and the Water Drop Penetration Time (WDPT) test [357]. These tests are based on measuring the infiltration time of a drop placed on the medium. More details on these tests are given in chapter 4. Results of these tests are related to contact angles and are more reproducible than direct measurements [188].

Regardless of the method used, we obtain wettability tendencies of the porous medium. Moreover, it is important to note that the wettability of a porous medium varies according to the amount of water in the medium [219]. As the water content decreases, the contact angle increases with the infiltration time in the porous medium. Therefore, a drier medium is more hydrophobic than a wet one.

Muehl et al. [241] define the notion of subcritical hydrophobicity, which occurs on perfectly dry surfaces, for $\Theta < \frac{\pi}{2}$.

However, in a dry medium, water is confined to a void volume and does not diffuse as a thin film around grains or aggregates, so the hydrophobic character is strongly accentuated. Beltrame [19] considers that hydrophobic interactions manifest at low saturation and when the medium saturates, it becomes hydrophilic again. This phenomenon can be observed using tomography equipment such as a confocal microscope [105] at pore scale, or drop infiltration measurements [162, 82, 73].

These angle variations are also impacted by the presence of amphiphilic matter. Chen et al. [77]
have shown that the contact angle increases with the concentration of amphiphilic matter for a water droplet placed on soil fractions lower than 2 µm. Moreover, Kroener et al. [201, 198] have revealed the existence of a critical mucilage concentration defining a percolation threshold. Thus, wettability varies in the presence of organic matter.

Whang and Wallach [354, 355] propose a model based on the time modification of the contact angle for a given concentration of amphiphilic matter. However, the models and measurements do not take into account the hydrophobic effect in the medium except by the contact angle. A porous medium must have hydrophobic interactions to generate partial wetting even if the medium is hydrophilic.

The question is how wettability is represented in the porous medium? Hydrophobicity and, more generally, the physico-chemical properties of the soil influence flow in macropores, and generally favor preferential flow [105, 19] (and pollutant transport). Then how can parameters linked to wettability be measured and integrated into soils? In order to provide some answers to this question in the chapter 2, we propose to review the wettability energy approach on a flat substrate and the models existing in the literature to model it.

### 1.3.3 Wettability energies

#### Disjoining pressure in plan surface

Disjoining pressure is a pressure representing the interactions between two interfaces that are separated by a distance $ζ$ and interact.

As presented in section 1.3.1, surface tension results from the forces acting at the molecular scale of an interface. For a volume of liquid flowing over a wall, the solid/liquid and liquid/gas interfaces form a surface tension associated with the liquid.

<figure><img src="image_22.png" /><figcaption>Figure 1.15: Sketch of the evolution of the free energy density of an infinitesimal infinitesimal film slice as a function of thickness (a) $γ_{film}(ζ) = γ_{SL} + γ$ and (b) $γ_{film}(ζ) = γ_{SG}$.</figcaption></figure>

However, in the neighbourhood of the contact line, it becomes more difficult to separate the two interfaces zones. The principle of disjoining pressure does not require an accurate description of the interface. It is therefore not necessary to model the effects of the triple line on the flow.

Consider two critical situations with a flat film of thickness $ζ$ on a wall:

i. the film is thick, i.e. $ζ$ is very large;

ii. the film is ultra-thin, i.e. $ζ$ is practically non-existent.

The energies in these two situations are different: in case (i), the film surface energy is $γ_{SL} + γ$ because the two interfaces zones are distinct from each other, whereas in case (ii), it is equivalent to the water
free surface tension $\gamma_{SG}$ (Figure 1.15). The surface energy associated with the water film must then be considered for a given film thickness.

To model this change in interface behavior, Derjaguin and Churaev [97] consider the two critical situations as without thickness interfaces and introduce an energy $V$ between the interfaces, which depends on the liquid thickness separating them. The surface energy associated with the flat water film can then be defined as a function of the film thickness:

$$ \gamma_{film}(\zeta) = \gamma_{SL} + \gamma + V(\zeta), \qquad (1.55) $$

where $V$ is a surface energy representing the continuity between the two situations. Thus, $V$ satisfies the following conditions:

$$ \bullet V(0) = \gamma_{SG} - \gamma_{SL} - \gamma = \mathcal{S}; $$

$$ \bullet \text{ interfaces do not interact and } \lim_{\zeta \to +\infty} V(\zeta) = 0. $$

Note that when the interfaces interact, i.e. the two zones partially overlap, a non-zero $V$ energy must be considered.

Derjaguin and Churaev [97] and De Gennes [92] then introduce the notion of "disjoining pressure" $\Pi(\zeta)$ associated with this energy. They considered disjoining pressure to be the negative of the derivative of the surface energy function $V$,

$$ \Pi = -\frac{dV}{d\zeta}. \qquad (1.56) $$

The surface energy $V$ is expressed as a function of the disjoining pressure for a thickness $\zeta$,

$$ V(\zeta) = \int_{\zeta}^{\infty} \Pi(e)de. \qquad (1.57) $$

The form of $V$ and $\Pi$ are used to express the wettability behavior between the two critical cases (i) and (ii). Different forms of function exist to represent the wetting of the drop on the solid [182].

<figure><img src="image_23.png" /><figcaption>Figure 1.16: Sketch of partial wetting (a) hydrophilic and (b) hydrophobic tendency. Red (resp. green) arrows represent repulsive (resp. attractive) interactions.</figcaption></figure>

At the mesoscopic scale between the solid surface and the free water surface, partial wetting occurs when repulsive interactions exist at least at short range. A competition between attractive (e.g. Van der Walls) and repulsive interactions can then be considered, as indicated by the green (resp. red) arrows representing the attractive (resp. repulsive) interactions in Figure 1.16). This representation
shows the two types of partial wetting: the hydrophilic tendency (Figure 1.16a) and the hydrophobic tendency (Figure 1.16b). So the shape of the surface energy V and the disjoining pressure Π influence the wettability of the liquid as a function of film thickness. Moreover, De Gennes [93] considers that if S > 0 (resp. S < 0) and V is decreasing (resp. increasing) at long range, then wetting is complete (resp. partial), as illustrated in Figure 1.17a (resp. Figure 1.17b).

<figure><img src="image_24.png" /><figcaption>Figure 1.17: Sketch of disjoining energy for (a) complete wetting, (b) partial wetting and (c) pseudo-partial wetting.</figcaption></figure>

Note that if $V(\zeta)$ decreases, the film will be thicker, and $\Pi(\zeta)$ is the pressure that must be applied to maintain the film at a thickness $\zeta$. Otherwise, $\Pi(\zeta)$ is the depression that must be imposed on the film surface to prevent it from vanishing the film thickness.

When the film thickness is zero, a singularity is observed at the triple line for partial wetting, as there is a discontinuity between an area of the dry solid substrate in contact with the drop. To avoid this problem, it is important to define a minimum film thickness that coexists<fnref n="2" /> with thicker films. Lallement [203] proposes to set a thickness $\zeta^*$ which avoids the singularity. However, this condition is too restrictive on the residual thickness. Other authors [93, 190] propose to define a precursor film $\zeta_p$ allowing to have a family of stable film<fnref n="3" /> ($\zeta \le \zeta_p$) coexist with thicker films. In order to avoid spontaneous dewetting of this equilibrium film, the disjoining pressure and its energy must be adapted such that:

$$ \begin{cases} V(\zeta_p) = \gamma_{SG} - \gamma_{SL} - \gamma, \\ \Pi(\zeta_p) = 0. \end{cases} \qquad (1.58) $$

This results in the increased Young-Dupré relation, relating the notion of contact angle and the disjoining energy:

$$ V(\zeta_p) = S = \gamma (\cos \Theta_e - 1). \qquad (1.59) $$

However, if there is not enough liquid to cover the entire surface for a film of thickness $\zeta_p$, a drop of liquid will spread out completely. In this case, two flat films of different thickness may coexist. In this case, we assume that the wetting is pseudo-partial rather than partial (Figure 1.18).

<footnote n="2">²Two films of thickness $\zeta_1$ and $\zeta_2$ coexist if and only if:</footnote>

<footnote n="1">1. $\Pi(\zeta_1) = \Pi(\zeta_2)$;</footnote>

<footnote n="2">2. $V(\zeta_2) = \Pi(\zeta_1)(\zeta_1 - \zeta_2) + V(\zeta_1)$.</footnote>

<footnote n="3">³A film of thickness $\zeta$ is stable if $P(\zeta)$ has a positive curvature ($P''(\zeta) > 0$)</footnote>
<figure><img src="image_25.png" /><figcaption>Figure 1.18: Sketch of pseudo-partial wetting.</figcaption></figure>

Therefore, the notion of disjoining pressure makes it possible to define different wetting models:

* complete wetting if $S > 0$ and $V$ decreasing at long-range;

* partial wetting if $S < 0$ and $V$ increasing at long-range;

* pseudo-partial wetting if $S > 0$ and $V$ increasing at long-range (Figure 1.17c).

## Pseudo-partial wetting model

Pseudo-partial wetting is wetting where the coexistence of a droplet and an ultra-thin precursor film occurs [328]. Experimentally, these precursor films have been observed in specific cases [163].

The precursor film is used as the *cut-off* parameter required to avoid the triple-line singularity of a wetting front [33]. Thus, the advancing and receding angles (angle derived from surface defects) do not defined with this model.

In recent years, this type of wetting model has been used to develop lubrication-type models based on the free energy of the system [271, 236, 257, 321]. These gradient-type models generally represent the dynamics of a drop on a substrate: bacterial biofilm [336, 337], liquid in the presence of surfactant [322] or on polymer brushes [323, 155, 165].

The question is how this type of modeling can be integrated into the soil? Cueto and Juanes [87, 86] by introducing a free energy per unit volume in the soil $\mathcal{E}$ as the sum of local potential energy (gravitational $\mathcal{E}_{gr}$ and capillary pressure $\mathcal{E}_{cap}$) and non-local potential energy $\mathcal{E}_{nl}$ [84]:

$$ \mathcal{E} = \mathcal{E}_{gr} + \mathcal{E}_{cap} + \mathcal{E}_{nl}. \qquad (1.60) $$

This energy purported by Cueto and Juanes does not take into account wettability notion and physico-chemical interaction. An interesting idea could be to relate this energy to the disjoining pressure with the saturation distribution in the soil. Thus, this notion of energy in the soil deserves to be explored in more detail to generalize water transfer dynamics, particularly in the presence of amphiphilic matter.

## 1.4 Conclusion

The main theoretical concepts relating to the effect of amphiphilic matter on the soil, capillarity and wetting are now defined. We have seen that the soil in the presence of amphiphilic matter includes the wetting mechanism using basic concepts such as surface tension and static contact angle. However, these concepts are insufficient to describe water transfer mechanisms, which become more complex in the presence of mucilage in the rhizosphere, for example.
The literature offers many possibilities for improving the representation of water transfers in soils in the presence of amphiphilic matter, but these remain open. The next two parts of this thesis are devoted to the modeling and numerical simulation of water transfers based on an energetic representation by analogy with the lubrication equation. An application will be realized to the rhizosphere in part III.
# Part II

## Model and numerical scheme
<figure><img src="image_26.png" /></figure>
# Chapter 2

## Wettability and free energy in soil

The aim of this chapter is to model water transfer in a soil in the presence of hydrophobic matter. We introduce a free energy associated to the transfer dynamics and we neglect inertial effects in the flow processes. This approach allows to build a model consistent with thermodynamics. One of the advances of this model is that the asymptotic dynamics are no longer unique, unlike current models. It is then possible to simulate hysteresis effects in soil flows. Finally, by modifying only the boundary conditions and the initial water content, we highlight water trapping phenomenon, fingering flows and intermittent-like dynamics.

This chapter is organised as follows. First, we describe the energy components in the soil in the presence of hydrophobic matter and define a free energy associated with these components in section 2.1. Water transfer modeling is presented in section 2.2 and its mathematical properties are described in sections 2.3-2.4. Finally, we present the water transfer simulation results from a published article (section 2.5).

## Contents

<table><tr><td><strong>2.1 Free energy in soil</strong></td><td><strong>40</strong></td></tr><tr><td>2.1.1 Thermodynamic approach</td><td>40</td></tr><tr><td>2.1.2 Free energy approach</td><td>42</td></tr><tr><td>2.1.3 Bulk energy</td><td>44</td></tr><tr><td>2.1.4 Wettability Energy</td><td>44</td></tr><tr><td>2.1.5 Water free surface energy</td><td>48</td></tr><tr><td>2.1.6 Conclusion</td><td>50</td></tr><tr><td><strong>2.2 Water transfer model</strong></td><td><strong>50</strong></td></tr><tr><td>2.2.1 Homogeneous porous medium</td><td>51</td></tr><tr><td>2.2.2 Stratified porous medium</td><td>53</td></tr><tr><td>2.2.3 Boundaries conditions</td><td>55</td></tr><tr><td>2.2.4 Dimensionless water transfer in homogeneous medium</td><td>56</td></tr><tr><td><strong>2.3 Theoretical results of homogeneous soil</strong></td><td><strong>57</strong></td></tr><tr><td>2.3.1 Wetting pressure</td><td>57</td></tr></table>
<table><tr><td>2.3.2</td><td>Equilibrium</td><td>58</td></tr><tr><td>2.3.3</td><td>Linear stability of the uniform steady-state</td><td>61</td></tr><tr><td><b>2.4</b></td><td><b>Instabilities in stratified media</b></td><td><b>64</b></td></tr><tr><td>2.4.1</td><td>Conclusion</td><td>65</td></tr><tr><td><b>2.5</b></td><td><b>Publication of the model</b></td><td><b>65</b></td></tr><tr><td>2.5.1</td><td>Letter journal: EPL</td><td>66</td></tr><tr><td>2.5.2</td><td>Summary of the article</td><td>66</td></tr><tr><td>2.5.3</td><td>Article</td><td>67</td></tr><tr><td><b>2.6</b></td><td><b>Conclusion</b></td><td><b>75</b></td></tr></table>

## 2.1 Free energy in soil

### 2.1.1 Thermodynamic approach

The concept of hydraulic potential in soil was first defined by Buckingham in 1907 [48]. He considered this potential to be the attractive energy present between soil and water: capillary potential. He suggests that this notion be used in conjunction with more general concepts such as thermodynamic potential or free energy. Later, Gardner [144] demonstrated how capillary potential depends on the volumetric water content $\theta$. This concept of capillary potential is used to represent the "soil/hydraulic potential". Water in the soil interacts with certain force fields that influence the potential of free water and bound water (water adsorbed to solid matrix). Thus, the hydraulic potential replaces the arbitrary classifications of water in the soil (gravity water, available water content, capillary water) introduced at the beginning of soil physics [170]. Therefore, an energetic classification was set up based on the potential energy of water in the soil, for instance gravity potential and capillary potential.

However, the term potential, as used in the literature on soils, has led to confusion because potential has the same physical dimension as pressure, which is not necessarily the thermodynamic potential. This confusion appears early with the experimental measurement which refers to the pressure and not the thermodynamic potential. A measuring system called the tensiometer was developed by Richards in 1931 [284] to measure this potential and define the retention curve used in Richards equation. However the retention curve refers to the pressure and not the thermodynamics potential. This confusion is source of ad hoc term (particularly in coupled system) because construction of the model is based on the pressure which is an intensive<fnref n="1" /> variable whereas it should be based on extensive<fnref n="2" /> variable as free energy. This leads to confusion between the $P_c$ capillary pressure and the derivative of the free energy functional, which have the same physical dimension.

Leverett [212] has shown that the negative of the capillary pressure is analogous to Lewis and Randall's "partial free energy" [213] and Gibbs "chemical potential" [223]. Thus, the capillary pressure is related to the free energy of the water per unit of pore space $F$. The free energy of the water per unit of pore space $f$ changes with the variation of water saturation $s$,

$$df = -P_c ds. \qquad (2.1)$$

<footnote n="1">¹An intensive property is a physical quantity whose value does not depend on the quantity of substance measured.</footnote>

<footnote n="2">²An extensive property is a physical quantity whose value is proportional to the size of the system it describes, or to the quantity of matter in the system.</footnote>
Therefore the difference in free energy of the water per unit of pore space between two saturation states $s_1$ and $s_2$ in the same volume is given by,

$$ \Delta f = f_2 - f_1 = - \int_{s_1}^{s_2} P_c ds. \qquad (2.2) $$

<figure><img src="image_27.png" /><figcaption>Figure 2.1: Panel (a) displays the volumes $V_t$, $V_f$ and $V_w$. $V_t$ is the elementary volume of soil. In this volume, the hashed region corresponds to the fluid: water or gas. $V_w$ is the water volume included in $V_w$. Panel (b) displays the variation $dV_w$ of $V_w$.</figcaption></figure>

Eq. (2.1) can be proven using a standard thermodynamic reasoning. Let us consider an elementary volume $V_t$ of soil. The pore space is noted $V_f$ which contains water volume $V_w$ (Figure 2.1a). We aim to express the change in free energy when the volume of water has varied by $dV_w$ (Figure 2.1b). We assume all transformations are isothermal. Since the number of water molecules $N$ varies, we consider the Landau free energy, also called grand potential $J$. We note, $\mathcal{U}$ the internal energy of water in the soil, $T$ the temperature $T$, $S$ the entropy and $\mu$ the chemical potential. We neglect the gas phase. By definition, the Landau free energy is

$$ J = \mathcal{U} - TS - \mu N. \qquad (2.3) $$

The variation of the Landau free energy during water exchange is

$$ dJ = d\mathcal{U} - SdT - TdS - \mu dN - Nd\mu. \qquad (2.4) $$

Since exchange of water $dN$ takes place at the same pressure and temperature, there is no variation in chemical potential: $d\mu = 0$. Under the isotherm assumption, the Landau free energy reads:

$$ dJ = d\mathcal{U} - TdS - \mu dN. \qquad (2.5) $$

From the first thermodynamics law,

$$ d\mathcal{U} = dW + dQ = -P_c dV_w + TdS + \mu dN. \qquad (2.6) $$

With the isotherm assumption,

$$ dJ = -P_c dV_w. \qquad (2.7) $$

As $V_w = sV_f$ then $dV_w = V_f ds$, then the free energy per unit of pore space is

$$ d(J/V_f) = -P_c ds. \qquad (2.8) $$
Considering the free energy of the whole soil assumed to be as a closed system (no water exchange), then the Landau free energy $J$ is the Helmholz free energy $F$. Then, by identifying $J/Vf$ to the Helmholtz free energy per unit of pore space $f$ in Eq. (2.8), we retrieve the relation Eq. (2.1).

In recent reviews, soil water potential is based on Gibbs free energy [369, 221, 372]. It contains different energy components depending on the processes studied (capillarity, gravity, ...). The terminology "total potential" is used to define this potential with the capillary pressure $P_c$. This definition is a source of confusion because the term potential refers to a thermodynamic state of a system, such as free energy. For instance, in Encyclopedia of Soils in the Environment, Tuller et al. [340] defines the total potential and its components as the thermodynamic potential. However, during the experimental application, this potential is identified as a measure of the capillary pressure!

In 1960, Philip [270] already indicates that the total potential $\Phi$ in literature corresponds to the Gibbs differential function. In our case, it will correspond to $-\frac{df}{ds} = P_c$. The author points out that it may be important to use the "absolute" rather than the differential form of thermodynamic function to study the evolution between two states. This absolute form is defined for a small element of soil with a given volume water content as equal to "minus the work required, isothermally and reversibly, per unit quantity of water, to extract all the water from the soil and transform it into pure, free water at a particular reference level" [269, 170]. Namely, the absolute form corresponds to the free energy variation $f_2 - f_1$ defined in Eq. 2.2.

### 2.1.2 Free energy approach

We propose to work using the notion of free energy, thus to avoid any ambiguity, we assume the free energy function $\mathcal{F}[s]$ on a domain $\Omega \subset \mathbb{R}^3$,

$$ \mathcal{F} = \int_{\Omega} f(s, \nabla s)d^3x, \qquad (2.9) $$

with the density free energy $f$ [J.m$^{-3}$] depending on $(s, \nabla s)$. In the rest of the manuscript, we will refer to the density free energy as the energy density.

The variation of the free energy are given by the Euler operator $L_f$ which is sometimes called the Lagrangian derivative of $f$. In physics, the expression $L_f(s)$ is called the variational derivative of the functional $\mathcal{F}$ at $s$. Thus the variation of the free energy are

$$ \delta\mathcal{F}(s, \varphi) = \int_{\Omega} \frac{\delta f}{\delta s} \delta s d^3x, \quad \forall \delta s := \varphi \in C_c^{\infty}(\Omega), \qquad (2.10) $$

where $C_c^{\infty}(\Omega)$ is the compact support of the domain $\Omega$.

This variational derivative of the functional $\mathcal{F}$ corresponds to distribution derivative,

$$ \delta\mathcal{F}(s, \varphi) = \langle \frac{\delta f}{\delta s}, \varphi \rangle. \qquad (2.11) $$

Here $\frac{\delta f}{\delta s}$ is a representative of the derivative of $f$ in the distribution sense. In general, we use the notation $\frac{\delta \mathcal{F}}{\delta s}$ instead of $L_f(s)$ or $\frac{\delta f}{\delta s}$. Note that the dimension of $\frac{\delta \mathcal{F}}{\delta s}$ is homogeneous with Joules per
unit volume and not Joules.

Using mathematical functional derivative, the variational derivative of the functional energy is related to the energy density,

$$ \frac{\delta \mathcal{F}}{\delta s} = \frac{\partial f}{\partial s} - \nabla \cdot \left( \frac{\partial f}{\partial \nabla s} \right). \qquad (2.12) $$

**Demonstration 2.1.1.** A functional derivative is derived from a non-normal directional differential [149]. It is defined by the relation

$$ \begin{aligned} \delta \mathcal{F}[s, \varphi] &= \left\langle \frac{\delta \mathcal{F}}{\delta s}, \varphi(s) \right\rangle, \\ &= \int_{\Omega} \frac{\delta \mathcal{F}}{\delta s} \varphi(s) dx, \\ &= \left[ \frac{d}{d\epsilon} \int_{\Omega} f(s + \epsilon\varphi, \nabla s + \epsilon\nabla\varphi) dx \right]_{\epsilon=0}, \qquad (2.13) \\ &= \int_{\Omega} \left[ \varphi \frac{\partial f}{\partial s} + \nabla\varphi \cdot \frac{\partial f}{\partial \nabla s} \right] dx, \\ &= \int_{\Omega} \varphi \left[ \frac{\partial f}{\partial s} - \nabla \cdot \left( \frac{\partial f}{\partial \nabla s} \right) \right] dx, \\ &= \left\langle \frac{\partial f}{\partial s} - \nabla \cdot \left( \frac{\partial f}{\partial \nabla s} \right), \varphi \right\rangle. \end{aligned} $$

Therefore,

$$ \frac{\delta \mathcal{F}}{\delta s} = \frac{\partial f}{\partial s} - \nabla \cdot \left( \frac{\partial f}{\partial \nabla s} \right). \qquad (2.14) $$

Thus the functional derivative of the free energy is the negative of the pressure in the system:

$$ \frac{\delta \mathcal{F}}{\delta s} = -P. \qquad (2.15) $$

It is therefore relevant to work on free energy and define the energy density associate to the different process in the soil then define the resulting pressure. This construction promotes the coupling of processes influencing water transfer and improve the modeled of the water transfer in the soil. When we have different processes, we can associate a free energy. Since the energy is an extensive variable, the global free energy is composed of different free energy. Each free energy are associated to an energy density. Therefore, we define the total energy density $f$ as the sum of the different free energy density,

$$ f = \sum f_{\square}, \qquad (2.16) $$

where the $\square$ symbolizes the different components of free energy.

**Example 2.1.2.** We consider two free energies $\mathcal{F}_1$ and $\mathcal{F}_2$ in a domain $\Omega$ such as the sum of the two free energy corresponds to the global free energy $\mathcal{F}$, i.e.

$$ \mathcal{F} = \mathcal{F}_1 + \mathcal{F}_2. \qquad (2.17) $$

We assume two density energies $f_1$ and $f_2$ associated to the free energy $\mathcal{F}_1$ and $\mathcal{F}_2$ respectively by the
relation Eq. (2.9). Thus the global free energy is written,

$$ 
\begin{aligned}
\mathcal{F} &= \int_{\Omega} f_1 d^3 x + \int_{\Omega} f_2 d^3 x, \\
&= \int_{\Omega} [f_1 + f_2] d^3 x, \\
&= \int_{\Omega} f d^3 x,
\end{aligned}
\qquad (2.18) $$

with the total energy density $f = f_1 + f_2$.

In the next sections, we describe the different contribution of free energy density $f_\square$. We distinguish bulk energy from wetting and surface energy.

### 2.1.3 Bulk energy

#### Gravity

Gravity is a fundamental force in physics, corresponding to the mutual attraction between two bodies of non-zero mass. In soil, this energy attracts liquids to the depths of the medium. It contributes to the downward flow. This term depends on the $z$ vertical coordinates. For an elementary $ds$ saturation variation, the elementary energy variation reads

$$ d f_{gr} = \rho g z ds. \qquad (2.19) $$

Note that for the study of horizontal transfer, this term varies only slightly and is therefore neglected in relation to variations in capillary energy.

#### Compressed air

In a porous medium close to the saturation, it is common to observe that air is trapped in the soil structure. In particular, the air is surrounded by water in micro-pores (Figure 2.2). The contact line may be pinned due to contact angle hysteresis. In this case the air is trapped in cavities of the soil structure. If this volume remains trapped, then variations in saturation result in a variation in air pressure. The negative of this compression work can be classically interpreted as free energy noted $f_{comp}$. Since this trapping phenomenon essentially only occurs for saturation close to one, this term depends on the saturation and is non-zero in a range of values close to saturation.

This term is also useful for numerical simulation, as it prevents the calculated saturation from exceeding unity. This allows us to successfully carry out numerical simulation of the phenomena of oversaturation at imbibition front and "infiltration refusals" in a porous layer.

### 2.1.4 Wettability Energy

Although the notion of wettability appears through various phenomena reported in the literature, such as the penetration time of a drop placed on a soil or the emergence of preferential flow, this concept is not explicitly defined in literature. According to Beltrame [19], we propose the following definition:
<figure><img src="image_28.png" /><figcaption>Figure 2.2: Sketch of air trapping in the soil at pore scale. Red dots represent pinned contact lines.</figcaption></figure>

**Definition 1. Soil wettability.**

*The wettability energy is the free energy due to the interaction between the porous matrix and water free surface.*

*Soil wettability refers to any property deriving from the wettability energy.*

In soil literature, this interaction usually refers to capillarity. Since Einstein's first paper in 1901 [117], it has been known that capillarity is due to interactions between a substrate and water. These interactions lead to a suction pressure and the capillary rise. Nevertheless, this only concerns attractive interactions and does not apply to hydrophobic matter.

Before we propose an analytical form of this energy, we review the physico-chemical origins of these interactions.

**Physico-chemical interactions**

In general, polar groups are hydrophilic. This is particularly the case with oxygen, nitrogen and sulfur groups. Hydrophilic properties are also influenced by hydrogen bonds. Thus, in the presence of OH groups, generally present in EPS, there are attractive interactions via these groups. In contrast, unpolar groups are not attracted to water. In particular, long carbon chains CH₃ and rings (e.g. aromatic).

<figure><img src="image_29.png" /><figcaption>Figure 2.3: Sketch of water molecule bounded to (a) a chloride ion through a hydrogen atom and (b) a sodium ion through an oxygen atom.</figcaption></figure>
EPS are long carbon chains with OH groups, so they are amphiphilic. Depending on their geometry in space, the molecule will be more or less hydrophilic.

Therefore attractive interaction depends on the composition and surface area of the porous matrix. In the case of water, its polar property is related to capillary action. This allows water to form hydrogen bonds using its partial charges. Water is a polar molecule due to the unequal sharing of electrons in the covalent bonds. Oxygen is more electro-negative and therefore brings the electrons closer to its nucleus. This creates a partial positive charge on the hydrogen atoms H and a partial negative charge on the oxygen atoms O. These partial charges allow water to form hydrogen bonds with other polar molecules. For example, water molecules will bind by their hydrogen atom to a particle containing chloride ions $Cl^-$ (anion negatively charged), while oxygen atoms will bind bond to a particle containing sodium ions $Na^+$ (cation positively charged) (Figure 2.3).

Therefore, water molecules adhere to the soil, they may be retained in the micro-pores of the soil. These interactions lead to attractive interaction in porous media. In contrast, the presence of organic matter can lead to repulsive interaction. Note that there are strongly hydrophobic interactions that form a coating on the surface of a solid, resulting in a non-wetting medium (e.g. TriMethylChloroSilane with sand or "magic sand" shown in Figure 2.4).

<figure><img src="image_30.png" /><figcaption>Figure 2.4: Photo of a drop on a sand with hydrophobic coating.</figcaption></figure>

## Modeling

In Brochard-Wyart et al. [41], substrate wettability is related to the different surfaces energies (see section 1.3). In a similar way, we can define a spreading parameter as the energy difference per aera unit $S = E_{dry}^m - E_{wet}^m$. Introducing the matrix surface tension $\gamma_P$ and the $\gamma_{PL}$ surface tension of a wetted porous matrix, then $S = (\gamma_P - \gamma_{PL})$. The soil is wettable if $S > 0$, it means that the energy decrease when it is wetted. In this case, capillary rise is observed. This phenomenon is common in porous media such as sugar in coffee or water rising in a sponge (Figure 2.5). The wettability energy compensates the gravity potential energy allowing a rise. It is customary in the literature to compare the porous medium to Jurin tubes. The height of capillary rise in soil is comparable to tubes of mean radius $r_m$ such as

$$h = \frac{2(\gamma_P - \gamma_{PL})}{\rho g r_m} \qquad (2.20)$$
<figure><img src="image_31.png" /><figcaption>Figure 2.5: Sketch of a capillary rise in a porous medium. The green arrow represents the rise of the liquid in the medium.</figcaption></figure>

Applying the Young equality: $(\gamma_P - \gamma_{PL}) = \gamma \cos \Theta$, where $\gamma$ is the surface tension of the liquid and $\Theta$ is static contact angle, the capillary rise then related to the contact static angle $\Theta$.

The drawback of this definition is that $\Theta$ is difficult to measure. Even if it could be measured for a substrate with the same physico-chemical properties, a Jurin tube does not behave in the same way as a porous medium. Some work has been carried out with corrugated tubes, but this is still far from reality in soils with very irregular, even fractal surfaces [267]. Finally, Jurin law does not take into account interactions between the free surface and the substrate. If the liquid film thickness $\zeta$ around the porous matrix is very thin, about a hundred nanometres then we have to take into account a $\zeta$-dependent surface energy. In low saturated porous media, these thin films certainly make a non-negligible contribution to the wettability energy. Then, for an element volume of wetted porous medium, we assume, in a first approximation, that the energy depends on the water content. The pressure in the porous medium derived from the energy density noted $f_w(s)$ is

**Definition 2.** *the wetting pressure defined as:*

$$ \Pi(s) = -\frac{df_w}{ds} $$

Assuming $f_w(1) = 0$, we get the integrated expression of $f_w$

$$ f_w(s) = \int_s^1 \Pi(s)ds $$

This function is analogous to the disjoining pressure for substrate wettability; we have used the $\Pi$ notation derived from Derjaguin. If we identify this pressure with the capillary pressure then it is a positive function and always decreasing. This implies that the energy $f_w(s)$ is decreasing ($\frac{df_w}{ds} < 0$) and convex ($\frac{d^2f_w}{ds^2} > 0$). In the context of liquid on a substrate, it is known that such a wetting free energy corresponds to total wetting, i.e. total spreading of the water. Mathematically, the convexity of the free energy implies the uniqueness of the solution, so the only equilibrium in a homogeneous medium is a homogeneous saturation distribution. This corresponds to Richards equation case, where the soil will moisten as much as possible to minimize its wettability energy. Note that this precludes
the existence of preferential flows as a steady state.

The presence of hydrophobic matter suggests that $\Pi(s)$ cannot always be decreasing, since water is repelled. We therefore propose to write $\Pi(s)$ as the sum of two functions:

$$ \Pi(s) = P_c(s) + P_r(s) \qquad (2.21) $$

the first term has the properties of capillary pressure: positive and decreasing. The second is negative and effectively represents repellent interactions. As for $P_c(s)$, the repellent pressure $P_r(s)$ vanishes near saturation. The wettability energy density writes

$$ f_w(s) = f_c(s) + f_r(s) $$

where

$$ f_c(s) = \int_s^1 P_c(s) ds \geq 0, $$

$$ f_r(s) = \int_s^1 P_r(s) ds \leq 0. $$

The wetting energy density is then the sum of the attractive energy part by the positive and decreasing $f_c(s)$ function and the repellent energy part by the negative $f_r$ function.

In very dry soils, residual water is observed experimentally. By analogy with pseudo-partial wetting with the existence of an ultra-thin film beyond the wetting line, we assume that for saturation close to zero, $\Pi$ is a decreasing function, i.e. $f_c \gg f_r$. In this case, we show in the next section that a residual saturation may remain, which we identify with the residual saturation discussed in the literature. For this reason, we have given the definition of saturation as the ratio of available volume to water volume and not the relationship $s = \frac{\theta_r - \theta_s}{\theta_s - \theta_r}$, where $\theta_r$ and $\theta_s$ are the residual and saturated water contents. There is therefore a possible shift in our curves when comparing retention curves. The advantage of our modeling is that $s_r$ and then $\theta_r$ are derived from the free energy and are not additional ad hoc parameters.

Note that we have introduced the repulsive term motivated by the presence of hydrophobic molecules. Nevertheless, the fact that the matrix is partially wetting, i.e. has an equilibrium contact angle $\Theta > 0$ implies that the disjoining pressure of a liquid on a flat substrate is necessarily increasing over a range of film thicknesses. Thus, in the case where very thin films predominate around the porous matrix, i.e. for low saturation, the disjoining pressure increases with film thickness. This indicates the existence of repulsive interaction whenever wetting is partial.

### 2.1.5 Water free surface energy

In the previous section, interaction energy between the free surface and the matrix is taken into account. However, the creation of a free surface also requires energy cost related to surface tension $\gamma$, regardless of interaction with the matrix. The energy cost to create the surface $A_f$ is $\gamma A_f$. Note, the air surface tension is negligible compared to the water surface tension. The energy density associated to this free surface is $f_\nabla = \frac{\gamma A_f}{\Omega}$, where $A_f$ is the water free surface in the elementary volume $\Omega$. In a continuous medium approach, it is difficult to estimate $f_\nabla$ and, above all, to relate this surface to
macroscopic quantities. The open is question is how to take into account this free energy?

A diffuse interface approach in multiphase media is proposed by Cahn Hilliard. In a first approximation, they showed that surface energy depends not only on concentration but also on the square of concentration gradient called "gradient energy" [54]. Instead to define a sharp interface between two fluids, a "phase field" $\rho$ varying continuously represents the presence of each fluid. For instance, when $\rho = 0$ the fluid is air while $\rho = 1$ the fluid is water. Any value in between represents the diffuse interface. The gradient of $\rho$ is large only in the diffuse interface of width $\epsilon$ (Figure 2.6).

It was proven that by minimizing the functional energy

$$ \mathcal{F}_{\epsilon} = \int_{\Omega} \frac{1}{\epsilon} f^{loc}(\rho) + \epsilon |\nabla \rho|^2 d\Omega, \qquad (2.22) $$

when $\epsilon \to 0$ (Gamma-convergence), the phase field $\rho(\mathbf{x})$ tends towards to the equilibrium distribution [237]. The minimal energy (inf($\mathcal{F}$)) is then the surface free energy. Recently, this functional is used as computation tool to track the interface of fluid mixing [8].

<figure><img src="image_32.png" /><figcaption>Figure 2.6: Surface tension in effective porous media.</figcaption></figure>

Cueto applied the Cahn-Hilliard approach to porous media where saturation is the equivalent of concentration in two-phase medium consisting of water and air. In other words, they make the assumption that the up-scaling at the elementary volume of porous medium has the same energy functional as the water/air interface. Then, the energy associated to the water/air interface expresses

$$ f_{\nabla} = \frac{1}{2} \Gamma |\nabla s|^2, \qquad (2.23) $$

with $\Gamma$ is the macroscopic surface tension coefficient. This physical quantity $\Gamma$ does not have the dimension of a surface tension (Energy per unit area) but its dimension is energy by unit length, i.e. a force in Newton. Note that in Cahn and Hilliard [54], the interface free energy depends also on the concentration $c$ what suggests that we should add a local term $f_{\nabla}^{loc}(s)$ in porous media. In the following, we associated the free surface energy only with the gradient term, called in Cueto and Juanes [86] non-local term since the local term can be included in the wetting free energy.

Adding the $f_{\nabla}$ energy to wetting free energy allows to raise the problem of ill-posed problem due to no convexity of $f_w$ as mentioned in Seppecher et al. [304]. The physical interpretation of mathematical regularizing role of $f_{\nabla}$ is discussed later.

Until Cueto-Juanes model, surface tension is not evoked in soil literature expected for wettability. But as we discussed in the beginning of the section the energy due to the interface is missing. This
expression is improved later in Beljadid et al. [17] to take into account the specificity of porous medium where the capillary energy interacts with the interface energy [17]:

$$f_{\nabla} \propto f_c(s) |\nabla s|^2. \qquad (2.24)$$

While the $|\nabla s|^2$ term proves crucial when instabilities are present (see Section 2.4), the improved expression Eq. (2.24) does not seem to bring new insights into the dynamics. Furthermore, the expression Eq. (2.24) for the derived pressure requires square roots of $f_c$, namely the wetting energy in our case. The latter can be negative in the case of highly hydrophobic media raising problem of its definition. Thus, in the present work, free surface terms will only be associated with the expression 2.23.

It is remarkable that none of the gradient-type expressions have been taken up in the literature outside the authors since its publication in 2008. A first explanation is that, apart from front instabilities (see Section 2.4), it does not provide fundamentally new behavior in the context of the Darcy-Richards equation. If the action of this gradient term is primarily to homogenize, as with capillarity, it has properties that cannot be described by simple capillarity. It prevents the formation of arbitrarily small or irregular structures. This is crucial during instabilities involving due to hydrophobic terms. Another, less obvious point is that it can generate flow front instability. It was the analogy with the viscosity instability of a falling film that motivated the introduction of this term by Cueto and Juanes [84].

### 2.1.6 Conclusion

To conclude, the density energy $f$ is the sum of four terms: wetting energy, compression energy, gravity energy and free surface energy. For convenience, in the rest of the manuscript, we include the compression energy term $f_{comp}$ in a "total wettability" term still noted $f_w$ and called wetting energy. With these notations, the density free energy reads

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

Regarding $f_w$, we distinguish the attractive from repellent interaction:

$$f_w = f_c + f_{rep} + f_{comp}, \qquad (2.26)$$

where $f_c$ is called capillary energy and it is associated to attractive interaction. The repellent part is noted $f_{rep}$. Note that if $f_{rep} \equiv 0$, then $\frac{df_w}{ds} = -P_c$ where $P_c$ has the same properties of capillary pressure presented in literature.

The expression of total wetting energy $f_w$ depends on the model and is specified in section 2.2.4.

## 2.2 Water transfer model

In this section, we present the development of a water transfer model based on the energy concepts mentioned in the previous section (Section 2.1). The mathematical model for a homogeneous medium is presented first, and then extended to heterogeneous media, in particular stratified media.
### 2.2.1 Homogeneous porous medium

In a homogeneous porous medium in the presence of hydrophobic matter, the water transfer in the medium $\Omega$ is modeled by describing its saturation profile $s(x, t)$ $[-]$ in a space $x \in \Omega$ during a finite time $t$ [s]. The dynamics associated to water movement is governed by the continuity equation (conservation law) of water density based on a hydraulic flow $J(s)$ [m.s$^{-1}$] (the s-conserving flux) such that,

$$ \phi \frac{\partial s}{\partial t} = -\nabla \cdot \mathcal{J}(s) = G(s), \qquad (2.27) $$

where the porosity $\phi$ $[-]$ is constant in the space during all the study time.

Now we need to relate the free energy $\mathcal{F}$ to the evolution equation of the saturation. We propose a gradient dynamics model which only requires the knowledge of the free energy, the flux of an order parameter (saturation in our case) and a quadratic dissipation functional. The method is based on the Onsager's variational principle [106, 107] as summarized in the appendix of [321].

In the following, we apply the method to the time evolution of the saturation. The total time derivative of the general energy functional $\mathcal{F}[s]$ is

$$ \frac{d\mathcal{F}}{dt} = \int_{\Omega} \frac{\delta\mathcal{F}}{\delta s} \frac{\partial s}{\partial t} d^3x, \qquad (2.28) $$

where $\Omega$ is the volume of soil studied and $\frac{\delta\mathcal{F}}{\delta s}$ is a derivative in the distribution sense. Injecting Eq. (2.27) in Eq. (2.28), the total time derivative of the general energy functional becomes,

$$ \frac{d\mathcal{F}}{dt} = - \int_{\Omega} \frac{\delta\mathcal{F}}{\delta s} \nabla \cdot \frac{\mathcal{J}(s)}{\phi} d^3x. \qquad (2.29) $$

We define the simplest dissipation functional such as quadratic in flux,

$$ \mathcal{W} = \int_{\Omega} \frac{1}{2} \mathcal{Z} \left| \frac{\mathcal{J}(s)}{\phi} \right|^2 d^3x, \qquad (2.30) $$

where the friction $\mathcal{Z} \ge 0$ may still depend on $s$. Using the Rayleighian [106, 321] $R = \frac{d\mathcal{F}}{dt} + \mathcal{W}$ with Eqs. (2.29)-(2.30), we obtain

$$ R = \int_{\Omega} \left[ \frac{1}{2} \mathcal{Z} \left| \frac{\mathcal{J}(s)}{\phi} \right|^2 - \frac{\delta \mathcal{F}}{\delta s} \nabla \cdot \frac{\mathcal{J}(s)}{\phi} \right] d^3x. \qquad (2.31) $$

By integrating by parts, Eq. (2.31) becomes,

$$ R = \int_{\Omega} \left[ \frac{1}{2} \mathcal{Z} \left| \frac{\mathcal{J}(s)}{\phi} \right|^2 + \frac{\mathcal{J}(s)}{\phi} \nabla \frac{\delta \mathcal{F}}{\delta s} \right] d^3x. \qquad (2.32) $$

We now minimize $R$, i.e we considered $\frac{\delta R}{\delta s} = 0$ and we obtain the flux,

$$ \mathcal{J}(s) = - \frac{\phi}{\mathcal{Z}} \nabla \frac{\delta \mathcal{F}}{\delta s}. \qquad (2.33) $$
The flow function $J(s)$ is proportional to the hydraulic conductivity in the medium $K(s)$ [Pa$^{-1}$.s$^{-1}$] and the energy functional $F[s]$,

$$J(s) = K(s)\nabla\frac{\delta F}{\delta s}, \qquad (2.34)$$

with conductivity function $K(s) = \frac{\phi}{\bar{z}} \ge 0$.

Therefore Eq. (2.27) can be rewrite,

$$\phi \frac{\partial s}{\partial t} = \nabla \cdot \left[ K(s) \nabla \frac{\delta F}{\delta s} \right]. \qquad (2.35)$$

In the soil, the hydraulic conductivity $K(s)$ is proportional to a conductivity at saturation $K_{sat}$ [m.s$^{-1}$] and a relative permeability $k_r$ [–],

$$K(s) = \frac{K_{sat}}{\rho g} k_r(s). \qquad (2.36)$$

By identification with the classical Darcy law, the functional derivative of the free energy $\frac{\delta F}{\delta s}$ is related to the pressure $P$ in the medium such that

$$\frac{\delta F}{\delta s} = -P. \qquad (2.37)$$

We have defined the energy density in the soil $f$, Eq. (2.25), using its components of gravity, wettability and surface free energies,

$$f(s, \nabla s) = f_w(s) + f_{\nabla}(s, \nabla s) + f_{gr}(s) = f_w(s) + \frac{\Gamma}{2} |\nabla s|^2 + f_{gr}(s). \qquad (2.38)$$

Then the functional derivative of the conservative part of the free energy $F[s]$ is

$$\frac{\delta F}{\delta s} = \frac{\partial f_w}{\partial s} - \Gamma \Delta s - \rho g z. \qquad (2.39)$$

Thus a gradient dynamics model based on the system governed by Eq. (2.27) is written,

$$\phi \partial_t s = \nabla \cdot \left[ K(s) \nabla \frac{\delta F}{\delta s} \right], \qquad (2.40a)$$

$$\frac{\delta F}{\delta s} = \partial_s f_w(s) - \Gamma \Delta s - \rho g z. \qquad (2.40b)$$

where $\Gamma$ is the so-called macroscopic surface tension [N]. Two fields are described in this system: the saturation evolution governed by Eq. (2.40a) and the pressure field Eq. (2.40b).

An analogy with models representing thin films on a substrate is given by this way of formulating the system governed by Eq. (2.40). In particular, it adopts the gradient dynamics formulation described in thin film literature [236, 257, 271, 321]. Thus, the dimensionless disjoining pressure $\Pi_w$ associated with partial wetting is defined as

$$\Pi_w(s) = - \frac{1}{\rho g H_{cap}} \partial_s f_w(s). \qquad (2.41)$$

Note that the dimensionless disjoining pressure is analogue to the J-Leverett function [212].
In this section, we define the dimensionless disjoining pressure $\Pi_w$ given by Eq. (2.41), using the pressure associated to the components of the wetting density energies such that

$$ \Pi_w(s) = h_c(s) + h_r(s) + h_{comp}(s), \qquad (2.42) $$

where $h_c$ is the dimensionless capillary pressure, $h_r$ is the dimensionless repellent pressure and $h_{comp}$ is the dimensionless compressed air pressure.

### 2.2.2 Stratified porous medium

A soil is a complex porous medium containing multiple heterogeneities. The thermodynamic model defined for a homogeneous medium, Eq. (2.40), is generalized to heterogeneous media. The dimensional formulation is a system of functions depending on space and time. In the case where the porosity is a constant in time then the system is given by:

$$ \left\{ \begin{aligned} \phi(x)\partial_t s &= \nabla \cdot \left[ K(s, \mathbf{x}) \nabla \frac{\delta F}{\delta s} \right], \\ \frac{\delta F}{\delta s} &= \partial_s f_w(s, \mathbf{x}) - \Gamma(\mathbf{x})\Delta s - \rho g z, \end{aligned} \right. \qquad (2.43) $$

where $\mathbf{x} = (x, y, z)$ are the spatial coordinates of the domain.

In the rest of this section, porosity is assumed to uniform porosity in each heterogeneous zone over time in each zone of the medium, i.e. $\phi(x) = \phi$.

A stratified medium is a medium made up of porous strata with different properties. In the most simple case, we consider two different strata in one dimension, where the transfer is vertical with $S_1 = [z_1, z_2]$ the domain of the first stratum and $S_2 =]z_2, z_3]$ the domain of the second stratum. The motion of the fluid in each stratum, according to Eq. (2.40), is governed by

**Stratum n°1:** for all $z \in S_1$,

$$ \phi(x)\partial_t s = \nabla \cdot \mathcal{J}_1 = \nabla \cdot [K_1(s)\nabla(\partial_s f_{w1}(s) - \Gamma_1\Delta s - z)], \qquad (2.44) $$

**Stratum n°2:** for all $z \in S_2$,

$$ \phi\partial_t s = \nabla \cdot \mathcal{J}_2 = \nabla \cdot [K_2(s)\nabla(\partial_s f_{w2}(s) - \Gamma_2\Delta s - z)], \qquad (2.45) $$

Parameters with a number index correspond to functions associated with the numbers of the stratum.

Eqs. (2.44)-(2.45), are naturally expressed as the same system represented by Eq. (2.43) using the formulation:

$$ \phi\partial_t s = \nabla \cdot [\mathcal{J}_1 \mathbf{1}_{[z_1, z_2]}(z)] + \nabla \cdot [\mathcal{J}_2 \mathbf{1}_{[z_2, z_3]}(z)], \qquad (2.46) $$

where $\mathbf{1}_X$ is the indicator function on the space $X$ such that

$$ \mathbf{1}_X(z) = \begin{cases} 1, & \text{if } z \in X, \\ 0, & \text{else.} \end{cases} \qquad (2.47) $$

However, it is not appropriate to use this indicator function because it introduces a discontinuity
between the two strata. A natural diffusion is possible between the two strata, in particular due to the mixing of particles between each stratum. There is a zone defined by the interval $[z_2 - \epsilon, z_2 + \epsilon]$, with a small $\epsilon$, where physically the two types of media coexist. A continuous heterogeneity function $f_{het}$ expresses the heterogeneity of the medium,

$$f_{het}(z) = \begin{cases} 1, & \text{if } z \in [z_1, z_2 - \epsilon], \\ 0, & \text{if } z \in [z_2 + \epsilon, z_3]. \end{cases} \qquad (2.48)$$

Figure 2.7 displays an example of a medium stratified using a function $f_{het}$.

<figure><img src="image_33.png" /><figcaption>Figure 2.7: Stratified porous medium composed by two strata.</figcaption></figure>

Stratum n°1 is considered as the reference stratum of the medium and stratum n°2 is expressed as a function of the characteristic functions of the reference stratum,

$$ \begin{aligned} K_2 &= (1 + C_K) K_1, \\ \Pi_2 &= (1 + C_{\Pi}) \Pi_1, \\ h_{r2} &= (1 + C_{rep}) h_{r1}, \\ \Gamma_2 &= (1 + C_{\Gamma}) \Gamma_1. \end{aligned} \qquad (2.49) $$

The $1 + C_{\square}$ are growth rates between the reference stratum and stratum n°2. The $\square$ index represents the function associated to the functions characterising the reference stratum.

Note that if all the constants $C_{\square}$ are zero then the homogeneous system is found.

We define various functions $f_{\square}$ as $f_{\square} = 1 + C_{\square}f_{het}(z)$ which allow us to express on the domain $S = S_1 \cup S_2$, a system similar to Eq. (2.43),

$$ \left\{ \begin{aligned} \phi \partial_t s &= \nabla \cdot \left[ f_K(z) K_1(s) \nabla \frac{\delta \mathcal{F}}{\delta s} \right], \\ \frac{\delta \mathcal{F}}{\delta s} &= \partial_s f_w(s, z) - f_{\Gamma}(z) \Gamma_1 \Delta s - \rho g z, \end{aligned} \right. \qquad (2.50) $$
with a wettability function

$$ \partial_s f_w(s, z) = -f_{\Pi}(z) \rho g H_{cap} \Pi_w(s, z) = -f_{\Pi}(z) \rho g H_{cap} [h_{c1}(s) + f_{rep}(z)h_{r1}(s) + h_{comp1}(s)]. \quad (2.51) $$

This heterogeneous system is studied in more detail, in section 2.5.3, following the analysis of the homogeneous system (section 2.3).

### 2.2.3 Boundaries conditions

The conditions on the boundary, called boundary conditions, are used to formulate the problem by imposing a physical situation corresponding to a phenomenon at the extremities of the domain. Thus system governed by Eq. (2.40) or Eq. (2.43) with an initial condition and boundary conditions forms a Cauchy problem.

Let $\Omega$ denote our domain and $\partial\Omega$ denote the boundaries of the domain. The boundaries of the domain are taken to be a continuous boundary and on this boundary there is a condition on each unknown, i.e. the saturation $s$ and the conservative energy functional $\frac{\delta F}{\delta s}$. There are four types of boundary conditions that we present for a given variable $u$: the Dirichlet condition, the Neumann condition, the periodic condition and the flux condition (Robin-Fourier condition).

Before to describe these conditions, the boundaries of the domain are considered to be expressed by the following partition:

$$ \partial\Omega = \Gamma_D \cup \Gamma_N \cup \Gamma_F \cup \Gamma_P. \quad (2.52) $$

Each partition corresponds to a boundary condition for $u$: $\Gamma_D$ the Dirichlet condition, $\Gamma_N$ the Neumann condition, $\Gamma_F$ the flux condition and $\Gamma_P$ the periodic condition.

Let $g$ a real function depending on space. For all positive $t$ and for all $\mathbf{x} \in \partial\Omega$,

1. **Dirichlet condition:** the value of $u$ is imposed for any $\mathbf{x} \in \Gamma_D$ such that

$$ u(t, \mathbf{x}) = g(\mathbf{x}). \quad (2.53) $$

Note that if $g \equiv 0$, the Dirichlet condition is said to be homogeneous.

2. **Neumann condition:** the first derivative is imposed for any $\mathbf{x} \in \Gamma_N$ such that

$$ \partial_n u(t, \mathbf{x}) = \frac{\partial u}{\partial \mathbf{n}}(t, \mathbf{x}) = g(\mathbf{x}), \quad (2.54) $$

where $\frac{\partial u}{\partial \mathbf{n}}(\mathbf{x}) = \nabla u(\mathbf{x}) \cdot \mathbf{n}$ with $\mathbf{n}$ the exterior normal vector. Note that if $g \equiv 0$, the Neumann condition is said to be homogeneous.

3. **Periodic condition:** the profile is continuous between two extremities of the domain. Let $\Gamma_1$ and $\Gamma_2$ two subsets of $\Gamma_P$ then for all $x_1 \in \Gamma_1$ and $x_2 \in \Gamma_2$,

$$ u(t, \mathbf{x}_1) = u(t, \mathbf{x}_2). \quad (2.55) $$
4. Flow condition: a flow $\mathcal{J}$ is imposed for any $\mathbf{x} \in \Gamma_F$ such that

$$ \mathcal{J}(\mathbf{x}) = K(u)\nabla H(u) = g(\mathbf{x}), \qquad (2.56) $$

where $K$ and $H$ are real functions of $u$.

In our case, there are two unknowns who may have different boundary conditions. There are different combinations of boundary conditions. In the 1D case, there are 17 combinations of possible boundary conditions (Appendix C). A wide range of physical situations can be studied. The choice of boundaries is important to represent an accurate physical phenomenon. Therefore, the compatibility of the boundary conditions with the problem being studied must be confirmed.

### 2.2.4 Dimensionless water transfer in homogeneous medium

In this section, we focus on the formulation of the homogeneous system. The dimensionless system is required to analyse the physical behavior of the system depending on its parameters. Using the dimensionless time scaling $T^* = \frac{\phi L^*}{K_{sat}}$ and the dimensionless space scaling $L^*$ such that the dimensional coordinates are related to the dimensionless coordinates

$$ \left\{ \begin{aligned} x &= L^* \tilde{x}, \\ y &= L^* \tilde{y}, \\ z &= L^* \tilde{z}, \\ t &= T^* \tilde{t}, \end{aligned} \right. \qquad (2.57) $$

the dimensionless system associated to Eq. (2.40) is given by

$$ \left\{ \begin{aligned} \partial_{\tilde{t}} s &= \tilde{\nabla} \cdot \left[ k_r(s) \tilde{\nabla} \frac{\delta \tilde{F}}{\delta s} \right], \\ \frac{\delta \tilde{F}}{\delta s} &= \tilde{H}_{cap} \partial_s \tilde{f}_w - T \tilde{\Delta} s - \tilde{z}. \end{aligned} \right. \qquad (2.58) $$

The dimensionless macroscopic surface tension is $T = \frac{\Gamma}{\rho g L^{*3}}$ and the dimensionless capillary height is $\tilde{H}_{cap} = \frac{H_{cap}}{L^*}$.

Note in the case where the energy density consider only the capillary energy, i.e the functional energy corresponds to $T = 0$ without $h_r$ and $h_{comp}$ define in Eq. (2.42), we retrieve the Darcy model. In the case where $T \neq 0$ and the energy density consider the capillary and the compressed air energies, we have the Cueto-Juanes model.

In the rest of this chapter, we work on the dimensionless problem. To simplify the notation, we remove the tilde symbol from the dimensionless parameters.

Moreover, we consider that the relative permeability is given by the Van Genuchten-Mualem model [343],

$$ k_r(s) = \sqrt{s} \left[ 1 - \left( 1 - s^{\frac{1}{m}} \right)^m \right]^2, \qquad (2.59) $$
the capillary pressure given by the Brooks and Corey model [42]

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

the repellent pressure given by a Gaussian function

$$h_r(s) = -\frac{1}{\sigma\sqrt{2\pi}H_{cap}} \exp\left[-\frac{(s-s_r)^2}{2\sigma^2}\right] \qquad (2.61)$$

and the compressed air pressure given by the Cueto and Juanes model [86], with $\kappa > 0$ a compressed coefficient,

$$h_{comp}(s) = sign(1 - \lambda) \exp[-\kappa(1 - s)] s^{-\frac{1}{\lambda}} \left(1 + \frac{\lambda}{\lambda - 1} \kappa s\right). \qquad (2.62)$$

The parameters $\lambda > 0$ and $0 \le m \le 1$ are shape parameters for the characteristic retention curve and the hydraulic conductivity. The parameters $s_r$ and $\sigma$ are repellent coefficients.

## 2.3 Theoretical results of homogeneous soil

In this section, theoretical results are presented about the modeling in homogeneous medium with hydrophobic matter. Particular attention is devoted to the definition of wetting in the soil, which varies according to the water content in the medium, and its influence on the transfer behavior and equilibria resulting from this modeling.

### 2.3.1 Wetting pressure

The notion of wetting pressure is generally not defined in the soil. We propose here to compare the structure of the wetting energy density and pressure associate to the equivalent wetting pressure of two main model of the literature: Darcy-Richards [284] and Cueto-Juanes [85] model.

In the Darcy-Richards model, the wetting pressure or the disjoining pressure $\Pi_w$ is represented by the capillary pressure $P_c$:

$$\Pi_w(s) = -P_c(s). \qquad (2.63)$$

In the Cueto-Juanes [85] model, the pressure is constructed from free energy associated with capillarity $f_c$ and air compression $f_{comp}$. These components of the free energy are related to capillary pressure and compressed air pressure $P_{comp}$, respectively. Then, we have the following equivalence with our model

$$\Pi_w(s) = P_c(s) + P_{comp}(s). \qquad (2.64)$$

Therefore, the Darcy-Richards and the Cueto-Juanes model generate wetting pressures associated with decreasing wetting energies (Figure 2.8a). The monotonic behavior of pressure implies the existence of a one-to-one relationship between pressure and saturation. These two wetting models only differ for saturation close to the saturation $s = 1$. The pressure defined by the Cueto-Juanes model decreases suddenly at negative pressures, whereas for the Darcy-Richards model the pressure remains positive. This difference defines an energy minimum close to saturation for the Cueto-Juanes model,
whereas in the Darcy the energy minimum occurs at $s = 1$. (Figure 2.8b).

<figure><img src="image_34.png" /><figcaption>Figure 2.8: (a) Pressure and (b) free wetting energy density of our model (solid blue line), Cueto and Juanes model (red solid line) and Darcy-Richards model (yellow dashed line). Parameters are $\lambda = 0.972$, $\sigma = 0.2$, $s_r = 0.05$ and $\kappa = 200$.</figcaption></figure>

In our modeling, we assume that the wetting energy density is the sum of attractive $f_c$ and repellent $f_{rep}$ density energies and the air compression $f_{comp}$ as previously with the Cueto-Juanes model. Thus our disjoining pressure is given by

$$ \Pi_w(s) = -(P_c(s) + P_{rep}(s) + P_{comp}(s)), \qquad (2.65) $$

with $P_{rep}$ the repellent pressure associated to the repellent energy density.

The wettability function in our model is, on the other hand, a non monotonous pressure since $P_{rep} \neq 0$ (Figure 2.8a). There is a range where the disjoining pressure increases. This increasing zone means that two extrema can be defined in addition to the unique extremum of other models (Figure 2.8b). This formulation of the wettability energy gives equilibrium states that the Darcy-Richards and Cueto-Juanes models cannot obtain for unsaturated media.

In the theoretical study of the model, we will keep the same disjoining function.

### 2.3.2 Equilibrium

A physical phenomenon is considered to be at equilibrium once it no longer changes in time. The phenomenon is a stationary state. In this study, saturation in the soil is at equilibrium if the time derivative of saturation is vanished, i.e.

$$ \partial_t s = 0. \qquad (2.66) $$

Based on Eq. (2.58) that means,

$$ \nabla \cdot \left[ k_r(s) \nabla \frac{\delta \mathcal{F}}{\delta s} \right] = 0. \qquad (2.67) $$

In a one-dimensional problem, there is a constant $C_0$ such that

$$ k_r(s) \nabla \left( \frac{\delta \mathcal{F}}{\delta s} - G_r z \right) = C_0, \qquad (2.68) $$
where $G_r$ is a parameter depending on the axis studied. For a vertical axis, $G_r = 1$, while for a horizontal axis, $G_r = 0$.

We study equilibria in the case where the domain is horizontal, i.e. $G_r = 0$. The equilibrium states is reached for the saturation profile in the soil are given by

$$k_r(s)\nabla\frac{\delta\mathcal{F}}{\delta s} = C_0. \qquad (2.69)$$

Assuming there is no flux boundary condition, then $C_0 = 0$. Since the hydraulic conductivity is a strictly positive function, then there is a constant $C_1$ such that

$$\frac{\delta\mathcal{F}}{\delta s} = C_1. \qquad (2.70)$$

This relation imply that

$$-\partial_s f_w + T\Delta s + C_1 = 0, \qquad (2.71)$$

where $C_1$ takes the role of a Lagrangian multiplier related to the mass conservation. The constant $C_1$ can also be related to the saturation $s_0$ such that $C_1 = \partial_s f_w(s_0) = -\Pi_w(s_0)$ since a uniform distribution $s \equiv s_0$ is a solution of Eq (2.71).

<figure><img src="image_35.png" /><figcaption>Figure 2.9: Disjoining pressure of our model. (a) Saturation tendency associated to the wetting energy. Gray area (on the right) corresponds to compressed air tendency, green area corresponds to hydrophilic tendency whereas red area corresponds to hydrophobic area. (b) Saturation equilibria depending on the Lagrange multiplier. Orange area corresponds to one saturation equilibrium, red area corresponds to three possible saturation equilibria and red dots corresponds to two possible saturation equilibria. Parameters are $\lambda = 0.972$, $\sigma = 0.2$, $s_r = 0.05$ and $\kappa = 200$.</figcaption></figure>

We then study the equilibria of uniform distributions with the disjoining pressure given in section 2.3.1. The disjoining pressure increases between $s_{cb} \approx 0.22$ and $s_{bu} \approx 0.58$. The increasing, resp. decreasing, range is highlighted in Figure 2.9a by the green, resp. red, area.

For a fixed constant $C_1$, we have three situations as shown in Figure 2.9b depending on the crossing number between the horizontal line $C_1$ and the $\Pi_w$ curve.

(i) one crossing: orange zone;

(ii) two crossing: dotted red,
(iii) three crossing: red zone).

In the situation (i), there is an unique equilibrium $s \equiv s_0$ possible associated to $\Pi_w(s_0) = -C_1$ in the ranges $]0, s_{cb}]$ or $[s_{cu}, 1]$. This solution is stable since $f_w''(s_0) > 0$ (Figure 2.10). In contrast, for the situations (ii) and (iii), there exist at least two equilibria whose can be stable or unstable.

<figure><img src="image_36.png" /><figcaption>Figure 2.10: Stability of uniform saturation profile. The profile is stable if $f_w(s) > 0$ and unstable if $f_w(s) < 0$. Parameters are $\lambda = 0.972$, $\sigma = 0.2$, $s_r = 0.05$ and $\kappa = 200$.</figcaption></figure>

The condition of coexistence of two different saturation is given by the equalities of pressure and of energy [41]:

$$ \bullet \ \Pi_w(s_{min}) = \Pi_w(s_{max}) = \Pi_0; $$

$$ \bullet \ s_{min}\Pi_w(s_{min}) + f_w(s_{min}) = s_{max}\Pi_w(s_{max}) + f_w(s_{max}). $$

These conditions are equivalent to consider that $s_{min}$ and $s_{max}$ have the same tangent given by

$$ f_w(s) = f_w(s_{min}) - (s - s_{min}) \Pi_0. \qquad (2.72) $$

<figure><img src="image_37.png" /><figcaption>Figure 2.11: Coexistence of saturation. Parameters are $\lambda = 0.972$, $\sigma = 0.2$, $s_r = 0.05$ and $\kappa = 200$.</figcaption></figure>

Figure 2.11 displays the coexistence of the saturation: $s_{min} \approx 0.13$ (almost dry soil) and $s_{max} \approx 0.95$ (almost saturated soil). Both equilibria are stable since near $s_{min}$ and $s_{max}$ the $f_w$ curve is convex.
Therefore, our wetting pressure allows the coexistence dry/wet soil at the equilibrium.

Note, there is a third uniform equilibrium possible $s \equiv s_u$ in the range $[s_{cb}, s_{cu}]$ which is unstable because $f_w(s_u)$ is concave. Linear instability means that an infinitesimal perturbation of the uniform profile leads to an exponential development of the amplitude of the perturbation. This phenomenon is similar to the dewetting phenomenon of thin film on a substrate. In our context, we expect coexistence of dry/wett soil regions. Assuming that in each dry and wet region the saturation is almost uniform then the saturation in dry, resp. wet, region correspond to $s_{min}$, resp. $s_{max}$, saturation discussed above.

In the next section, we explain how to determine linear instability for equilibrium.

### 2.3.3 Linear stability of the uniform steady-state

In a homogeneous medium, due to the repulsive term $h_r$, the uniform distribution of water $s \equiv s_0$ can be linearly unstable with respect to infinitesimal perturbation. The stability analysis was already carried out in the lubrication framework of a thin film, for instance Thiele [190], or in a similar framework in Beltrame chapter [19].

Let $s_0$ a fixed saturation of Eq. (2.70) such that $s_0$ satisfy $\dot{G}(s_0) = 0$ where G is defined in Eq. (2.27). The linear stability of this dynamic is given using the linearize expression of Eq. (2.58) with $\epsilon \ll 1$ yields. Let a small $\epsilon$ and $u$ an integrable function in $\Omega$ such that $s(x) = s_0 + \epsilon u$. Inserting this expression in Eq. (2.58), for a fixed saturation $s_0$,

$$ \epsilon \partial_t u = -\nabla \cdot (K(s_0 + \epsilon u) \nabla (\Pi_w(s_0 + \epsilon u) + T \Delta(s_0 + \epsilon u))). \quad (2.73) $$

We linearize K and $\Pi_w$ thanks to Taylor expansion applied to $s_0$ with an error order of $O(\epsilon^2 u^2)$.

$$ K(s_0 + \epsilon u) = K(s_0) + \epsilon u \mathbf{D}_s K(s_0) + \text{h.o.t.} \quad (2.74) $$

and

$$ \Pi_w(s_0 + \epsilon u) = \Pi_w(s_0) + \epsilon u \mathbf{D}_s \Pi_w(s_0) + \text{h.o.t.}, \quad (2.75) $$

where D is the differential operator. Thanks to Eqs. (2.74)-(2.75) and the fact that we have Eqs. (2.70)-(2.73) becomes,

$$ \partial_t u = -\nabla \cdot [(K(s_0) + \epsilon u D_s K(s_0)) \nabla (u \mathbf{D}_s \Pi_w(s_0)) + T \Delta u)]. $$

Thus, when $\epsilon$ goes to zero we have:

$$ \partial_t u = -\nabla \cdot [K(s_0) \nabla (u \mathbf{D}_s \Pi_w(s_0)) + T \Delta u]. \quad (2.76) $$

The eigenvalues of the linear operator on the right hand determines the linear stability. If all the eigenvalues are negative the steady-state is stable otherwise it is unstable. The determination of these eigenvalues requires numerical tools presented in the next chapter 3 except if $s_0$ is a uniform distribution. We focus, in the following, to this particular case to analyse the stability.
The linear dynamics Eq (2.76) becomes for a uniform $s_0$ steady-state

$$ \partial_t u = -K(s_0)\Delta(u\mathbf{D}_s\Pi_w(s_0)) + T\Delta u. \qquad (2.77) $$

To assess the linear stability of the flat flow $s(x) = s_0$, we use a Fourier mode decomposition from $s(x) = s_0 + \epsilon \exp(\beta t + ikx)$ where $\beta$ and $k$ provide the growth rate (own eigenvalue) and the wave number of disturbances. Concretely, we take $u(X, t) = \exp(\beta t + ikX)$ in Eq. (2.77) and we derive the dispersion relation

$$ \beta(k) = K(s_0)(\mathbf{D}_s\Pi_w(s_0)) - Tk^2)k^2. \qquad (2.78) $$

We called $k_{crit}$ the critical eigenvalue such that

$$ \beta_{crit} = 0, $$

and if $T \neq 0$,

$$ k_{crit} = \sqrt{\frac{\mathbf{D}_s\Pi_w(s_0)}{T}}. \qquad (2.79) $$

If $0 < k < k_{crit}$ then $s_0$ is unstable whereas if $k > k_{crit}$ then $s_0$ is stable, since the linear stability of $s_0$ depends of the sign of $\beta$.

<figure><img src="image_38.png" /><figcaption>Figure 2.12: Critical mode depending on the saturation for our model, Cueto and Juanes model and Darcy-Richards model. Parameters are same parameters than Figure 2.8 with $T = 100$, $n = 3.31$ and $m = 1 - \frac{1}{n}$.</figcaption></figure>

Figure 2.12 displays the maximal mode $k_{crit}$ associated to our model, Cueto-Juanes model and Darcy-Richards model depending on the uniform saturation $s$. In the Cueto-Juanes and Darcy-Richards models the critical mode $k_{crit}$ is always zero. Indeed, this mode is due to conservative form of the governing equation and so all other modes $k$ are associated to a negative real part of the $\beta(k)$ eigenvalue. Thus, they do not correspond to a linear instability. On the contrary, our model generates instability in the saturation range $[s_{cb}, s_{cu}]$. The critical mode $k_{crit}$ rate reaches a maximum $k_M = 0.32$ for $s_{crit} = 0.31$. The maximum value $K_M$ indicates that there is minimum wavelength $\lambda_M = \frac{2\pi}{k_M}$ of unstable perturbations.
Now, we focus on the dispersion relation given the growth rate $\beta(k)$ w.r.t $k$ for our model. Figure 2.13 displays dispersion relation for different $s$ saturation values according to different surface tension coefficients $T$. Modes are stable if the eigenvalues $\beta$ associated to these modes are negative. In the cases $s = 0.1$ (Figure 2.13a) and $s = 0.7$ (Figure 2.13d), there are no critical modes then all the modes $k$ are stable. In the cases $s = 0.3$ (Figure 2.13b) and $s = 0.5$ (Figure 2.13c), there are two cases depending if $T = 0$ or if $T > 0$.

For $T > 0$, there exists a critical mode such that the eigenvalues $\beta$ associated to the modes between $k = 0$ and $k = k_{crit}$ are positive, otherwise the eigenvalues are negative.

For $T = 0$, all eigenvalues associated to the modes of the system are positive. Moreover the growth rate is not bounded when $k$ tends towards infinity. This means that infinitely small structures develop, which is not consistent with a macroscopic approach. The stability analysis points out that a finite tension surface coefficient is required to get a well-defined model by giving a cut-off scale of the unstable modes. This conclusion is consistent with theoretical results when a $f_w$ energy is non convex [304]. In the following, we assume that $T > 0$.

<figure><img src="image_39.png" /><figcaption>Figure 2.13: Dispersion relation of our model according to different surface tension constant $T$.</figcaption></figure>

The unstable modes are related to the Fourier modes. Each $k$ Fourier mode is associated to the wavelength $\lambda^*$ such that

$$k = \frac{2\pi}{\lambda^*}.$$

That means, there exists a critical wavelength $\lambda_{crit}^*$ which is related to surface tension and disjoining
pressure for which the uniform distribution is unstable

$$ \lambda_{crit}^* = 2\pi \sqrt{\frac{T}{\mathbf{D}_s \Pi_w(s_0)}} $$

All perturbations with a wavelength above the critical value $\lambda_c^*$ will develop exponentially. The resulting water distribution is a succession of dry and wet regions, the typical $\lambda_m^*$ size of which is related to the mode of fastest growth rate. According to Figures 2.13b and 2.13c there exists an unique maximum growth rate $\beta(k)$. Indeed, we can prove that the maximum $\beta(k)$ given by Eq. (2.78) is related to $\lambda_c^*$ by the relation

$$ \lambda_m^* = \sqrt{2\lambda_c^*} $$

In conclusion, if the mean saturation is in the hydrophilic tendency then the equilibrium profile is stable. Otherwise, the uniform profile is unstable only if the domain of study is larger than $\lambda^*$. In this case, an inhomogeneous saturation distribution is expected. This equilibrium is similar to the pseudo-partial wetting where an ultra-thin precursor film coexists with drops or thicker films. Here, the low saturation in dry region plays the role of the precursor film.

Note the importance of boundary conditions, since for domains smaller than $\lambda^*$ the uniform distribution is stable.

## 2.4 Instabilities in stratified media

Let us consider a vertical stratification of hydrophilic and hydrophobic layers as presented in section 2.2.2. According to experimental observation, water first infiltrates in the hydrophilic layer and then either stay trapped or for enough water content, or the water infiltrates into the hydrophobic layer via fingered [105, 307]. Since work in the seventies of [262] devoted to the study of the imbibition front instability, this problem has been the subject of many studies in order to understand and model this phenomenon in homogeneous or heterogeneous media. However, this problem still remains an open question as shown by Fürst et al. [140]. It is known that the Richards equation is unable to reproduce this instability even by varying the capillary height in each stratum [100]. Indeed the asymptotic dynamics is always a stationary flow where the water saturation is modulated by the heterogeneity. Many variants of the Richards equation, such as dynamic capillary pressure [168] and hypo-diffusion [119] have been developed to obtain out-of-equilibrium dynamics. However, the quantitative results are not in agreement with experiment, especially concerning the oversaturation observed at the infiltration front. Cueto-Felgueroso and Juanes [86] introduced the free surface energy discussed in section 2.1.5. The front instability then appears as a Saffman-Taylor type instability in homogeneous medium. Overshoot of the saturation or oversaturation at the front appears as a consequence of the energy term $|\nabla s|^2$. The oversaturation induces a greater mobility, leading to fingering in a homogeneous porous medium. However, this model does not capture the instabilities of an initially blocked front, for instance, at the interface of a hydrophilic layer on a hydrophobic layer.

One of the hypotheses of the link between hydrophobicity and fingering is the existence of the hysteresis of the contact angle which is revealed through the hysteresis of the infiltration curves [168, 16]. Rätz and Schweizer [276] showed the existence of an unstable front due to a contact angle hysteresis
introduced directly into a Richards-type equation. Concretely, in order to advance the "charge" must pass a threshold such as a depinning phenomenon.

In our model, at the imbibition front the repellent term are dominant and then a similar charge must pass a threshold arriving in the hydrophobic medium. Indeed, the governing equation of saturation is similar to the lubrication equation of the liquid bridge height driven by gravity on heterogeneous substrate stripped by periodic hydrophilic and hydrophobic coating. This configuration was deeply studied by Thiele and Knobloch [325, 324] then by Beltrame et al. [22] in 2D and 3D geometry. The case of liquid ridge blocked by an hydrophobic defect displayed very complex scenarios [21]. Under appropriate conditions the ridges may either depin and slide downslope as the slope is increased, or first break up into drops via a transverse instability, prior to depinning. The first instability is associated with the saddle-node bifurcation and is referred to in this context as the depinning instability. The resulting dynamics is an intermittent-like behavior, alternating a long period of blockage by the hydrophobic defect and relatively rapid flow after depinning. The second one is associated to a pitchfork bifurcation and it is called the Rayleigh-Plateau instability. This instability corresponds to the development of front undulations.

Both behaviors are reproduced by time integration and bifurcation analysis in the case of a stratified soil and presented in our publication [20]. At low water content, water is trapped in the hydrophilic layer. Above a critical water content, intermittent flow alternates between a long duration in the hydrophilic layer and rapid flow in the hydrophobic stratum. Prior this threshold, a Plateau-Rayleigh instability occurs, leading to the emergence of fingers, resulting from asymptotically stationary flow in the hydrophobic layer. The development of water distribution at the hydrophilic/hydrophobic interface is strongly reminiscent of the experiment description in Sililo and Tellam [307].

### 2.4.1 Conclusion

This modeling based on a free energy is suitable for studying problems in both homogeneous and heterogeneous media. The model allows several equilibria to coexist, in contrast to the single state with the water transfer approaches presented in literature. Indeed some of these states of equilibrium cannot be reproduced by other models. It is difficult to follow and predict the dynamics exactly. Indeed, there is no usual analytical expression due to the non-linearity of the system. In this situation, it is important to approximate the time evolution of saturation using numerical simulation. More details about the numerical tools used to simulate the problem are presented in the chapter 3.

## 2.5 Publication of the model

The work on the model and the simulation results presented in this chapter represent innovative work in the understanding of water transfer in a porous medium, especially in the presence of hydrophobic matter. These results were published in a letter of 7 pages (limit imposed by the journal) in the EPL-Europhysics Letters, in June 2022.

Firstly, the publication strategy and the choice of this letter journal are discussed in section 2.5.1. A summary of the main results of this modeling in a homogeneous porous medium and in a stratified medium are presented in this article with some additional simulations in section 2.5.2. The preprint
version of the article is introduced in section 2.5.3 and finally a global discussion of this work applied
to a soil is presented in section 2.6.

**2.5.1 Letter journal: EPL**

The article has been accepted and published in the letter journal EPL, edited by IOPscience. This internationally peer-reviewed scientific journal publishes letters containing work of great interest and importance to one or more sections of the physics community, including *Fluid and nonlinear dynamics*. The main purpose of publishing in this journal is to introduce this water transfer model into the physics literature, particularly in the field of porous media with nonlinear dynamics. Once this publication is in the literature, it will be easier to develop future work in more specific and complex porous media such as soils. In particular, the variety of transfer dynamics reproduced by the model can be compared qualitatively with experiments in the literature.

**2.5.2 Summary of the article**

In the literature and in experiments in porous media, various water transfers occur, such as water trapping, fingering and uniform transfers when the medium is sufficiently saturated. Some of these phenomena are reproduced in mathematical models, as presented in section 2.2. However, current models are unable to reproduce all these phenomena in a same model.

We return to the design of a wettability model based on the concept of free energy. This model is applied to a homogeneous medium and then to a stratified sandstone medium in the presence of a hydrophobic matter.

Numerical simulations revealed the formation of partial wetting in the homogeneous medium, with a profile that varies according to the mean saturation. In addition, a potential hysteresis dynamic is highlighted by a bifurcation diagram (see section 3.5) whose parameter is the mean saturation of the water in the medium. These first results show the existence of multiple equilibrium states varying according to saturation between partial wetting and total wetting.

In the stratified medium consisting of a hydrophilic stratum and a hydrophobic stratum repeated periodically, we revealed different transfer in 1D and 2D by varying only the mean saturation in the medium. In 1D, the numerical simulations captured:

* the dynamics of water trapping in the hydrophilic stratum at low saturation levels;

* the dynamics of intermittence in the hydrophobic stratum for saturation in a saturation hydrophobic range (see section 2.3.1);

* the uniform transfer dynamics produced by the Darcy-Richard model with high saturation.

These results are confirmed in 2D where, in addition to the intermittence phenomenon, we have the formation of fingering in the hydrophobic stratum with a blocked wettability front in the hydrophilic stratum. This fingering is in agreement with the experiments of Nieber et al. [250] and Bauters et al. [15]. Figure 2.14 represents the different dynamics obtained by 2D simulations of the water transfer model in a porous medium.
<figure><img src="image_40.png" /><figcaption>Figure 2.14: 2D water transfers in a stratified porous medium with a hydrophilic stratum on the top and a hydrophobic stratum on the bottom and different mean saturation in the medium.</figcaption></figure>

This model is able to reproduce several water transfer phenomena by only changing the mean saturation of the water in the porous medium. These results are promising for a better understanding of water transfer in the presence of hydrophobic substances and for taking into account the hysteresis phenomena acting on the water in the medium.

### 2.5.3 Article

The version of the paper proposed in this manuscript is the version submitted to the journal EPL-Europhysics Letters, Volume 138, Number 5 [20]. The published version can be found on the IOPscience website. The preprint version of the article appears on the following pages.
# Model of hydrophobic porous media applied to stratified media: Water trapping, intermittent flow and fingering instability

P. BELTRAME<fnref n="1" /> and F. CAJOT<fnref n="1" />,<fnref n="2" />

**Abstract** – A macroscopic model of water infiltration is developed and applied to a stratified porous medium containing hydrophilic and hydrophobic layers. The governing equation is a gradient dynamics as in the Richards model but in which a wettability term replaces the suction pressure, taking into account the effective action of the hydrophobic matter. The resulting governing equation is a fourth order PDE similar to the lubrication model or Cahn-Hilliard equations known to reproduce complex dynamics such as intermittency and fingering instability. The model makes it possible to simulate the finite spreading of water in horizontal directions when the porous medium is hydrophobic. Applied to a stratified medium alternating hydrophilic and hydrophobic layers, we highlight, depending on the mean water content, phenomena of water trapping in hydrophilic layers, intermittent flow and fingering instability. These regimes are studied using time integration bifurcation and numerical continuation tools.

**Introduction**. – Water transfer in unsaturated porous media displays phenomena that are still poorly understood such as fingering instability of the wetting front, water trapping and intermittent flow in heterogeneous media. These phenomena are enhanced in the presence of hydrophobicity due, for instance, to organic matter. Experimental observation shows that the water infiltrates into the hydrophobic layer via fingered flows even if the layer is homogeneous [1]. Then, the fingered flows result from an imbibition front instability. In Darcy-Richards modelling, the wetting property is defined via capillary rise and the pressure head. In such a model, the water content diffuses towards the drier region and thus can reproduce neither the water trapping nor the fingering instability [2].

Since work in the seventies of [3,4] devoted to the study of the imbibition front instability, this problem has been the subject of many studies in order to understand and model this phenomenon in homogeneous or heterogeneous media. However, this problem still remains an open question as shown by [2]. It is known that the Richards equation is unable to reproduce this instability. Many variants of the Richards equation, such as dynamic capillary pressure [5] and hypo-diffusion [6] have been developed to obtain out-of-equilibrium dynamics. However, the quantitative results are not in agreement with experiment, especially concerning the oversaturation observed at the infiltration front. Cueto-Felgueroso and Juanes [7] introduced

a surface energy in unsaturated media. The front instability then appears as a Saffman-Taylor type instability. Oversaturation at the front induces a greater mobility, leading to fingering in a homogeneous porous medium. However, this model does not capture the instabilities of an initially blocked front, for instance, at the interface of a hydrophilic layer on a hydrophobic layer. To take into account the front hysteresis, models based on the analogy of the contact angle of a water front on the surface are found in the literature, for instance in [8].

The hydrophobicity of the soil is another factor inducing instabilities which is not taken into account by classic models [9–11]. One of the hypotheses of the link between hydrophobicity and digitation is the existence of the hysteresis of the contact angle which is revealed through the hysteresis of the infiltration curves [5, 12]. Rätz and Schweizer [8] showed the existence of an unstable front due to a contact angle hysteresis introduced directly into a Richards-type equation. Concretely, in order to advance the "charge" must pass a threshold such as a depinning phenomenon. The role of contact angle hysteresis and hydrophobicity in front instability phenomena led to focus on wettability in porous media. In contrast, in the present study, we aim at describing the wettability via surface energy at the interfaces. An effective term of short-range repulsion is introduced to model the interaction between the porous matrix and thin water films. According to [13],

<footnote n="1">¹ Avignon Université – UMR1114 EMMAH INRAE-AU - F84000 Avignon - FRANCE</footnote>

<footnote n="2">² INRAE - UMR1114 EMMAH INRAE-AU - F84000 Avignon - FRANCE</footnote>
<figure><img src="image_41.png" /><figcaption>Fig. 1: Suction profile versus saturation. The fine line is the Genuchten-Mualem suction with $n = 10$ (Eq. 3), the bold line is the pressure $\Pi$ taking into account the hydrophobicity given by Eqs (4-5). Parameters are $h_f = 0.9$, $\delta = 0.1$, $s_{rg} = 0.2$.</figcaption></figure>

this term is non negligible at low saturation. The sum of this repulsive term with the capillary pressure constitutes the analogous of con- and dis-joining pressures describing the interaction of a partial wetting liquid [14]. An effective surface tension proposed in [7, 15] is also taken into account. After detailing the model, we apply the model to stratified media with different wettabilities.

Modelling. –

*Governing equation.* Let us consider a homogeneous porous medium with porosity $\phi$, capillary height $H_c$ and the conductivity $K_s$ at saturation. If the medium is unsaturated, we define the saturation $s \in [0; 1]$ as the ratio between the water volume over the interpore volume for a dry medium. The dimensionless Darcy-Richards governing equation reads:

$$ \frac{\partial s}{\partial t} = -\nabla \cdot \vec{\Phi}(s), \qquad (1) $$

where the flux $\vec{\Phi}$ is proportional to the relative permeability $k_r(s)$ and the pressure gradient due to the static pressure and the capillary pressure $h_c(s)$:

$$ \vec{\Phi} = k_r(s) \nabla (h_c(s) + z) \qquad (2) $$

with positive z-axis downwards. The space is scaled by $H_c$, the time by $\phi H_c/K_s$ and the pressure by $\rho g H_c$, $\rho$ is water density. In a hydrophilic medium, the dimensionless suction $h_c(s)$ is a positive and decreasing function. In this paper, we use the Genuchten Mualem model [16]

$$ h_c(s) = (s^{\frac{n-1}{n}} - 1)^{1/n}, \qquad (3) $$

with n a parameter taking the values typically in the range [1, 10] depending on the soil type. In this paper, we fixed the value to $n = 10$ which is relevant for a hygiene sandstone and widely used in the literature [15, 17]. Because $h_c(s)$ suction is a decreasing function of the saturation, the water diffuses from high saturation to low

saturation in a horizontal thin layer. Then, the final state of water imbibition is a uniform distribution $s(x) \equiv s_m$, where $s_m$ is the mean water content. In the presence of hydrophobic matter, water diffusion may be stopped due to the hydrophobicity, thus, the Darcy-Richards equation does not apply. The latter behavior is analogous to partial wetting of a fluid on a substrate. The droplet does not spread totally and makes a finite contact angle. The contact angle results from competition between attractive and repulsive interaction acting at different scales between the substrate and the free surface. Usually, the repulsive term acts at short range and the attractive term at longer range [18, 19]. These interactions define the substrate wettability. Our approach is to introduce a similar wettability term in porous media.

For unsaturated porous media, there is an interaction between the free surface of water and the matrix that is represented by the suction $h_c(s)$. If hydrophobic matter is present, we modify the potential $h_c(s)$ by $h_c(s) - h_r(s)$ where the repellent term $h_r > 0$ represents the repellent action. This term depends on the concentration of hydrophobic molecules but also on water saturation [13]. Indeed, if the hydrophobic molecules are surrounded by water, their hydrophobicity is weaker. Then, the repellent term $h_r(s)$ should be a strongly decreasing function. We chose an exponential function:

$$ h_r(s) = h_f \exp(-\frac{s}{s_{rg}}) + \delta, \qquad (4) $$

where $h_f$ is the magnitude of hydrophobicity, $s_{rg}$ is the saturation range affected by the hydrophobicity and $\delta > 0$ stands for the remaining repellent pressure at saturation. The resulting dimensionless pressure is then

$$ \Pi(s) = h_c(s) - h_r(s). \qquad (5) $$

In this paper, we fixed $s_{rg} = 0.2$, $h_f = 0.9$, $\delta = 0.1$, for hydrophobic media. The delta value means that the suction is reduced by 10% close to saturation. The values $h_f$ and $s_{rg}$ were chosen in order to induce fingering instability phenomena. The resulting curve $\Pi(s)$ is shown in Fig. 1. The function is no longer monotonic: $\Pi(s)$ increases on the $[s_1; s_2] = [0.045; 0.471]$ interval (see Fig. 1). Then, in the range bounded by the extrema $[\Pi(s_1); \Pi(s_2)] = [0.588; 0.842]$ there is no longer uniqueness of the saturation for a given suction $\Pi_0$. The consequence is that an equilibrium exists between low saturation which has a value in the interval $[0.0045; 0.0449]$ and high saturation in the interval $[0.471, 0.977]$ allowing a finite water diffusion in contrast to Darcy-Richards modeling. As a result a uniform profile $s = s_m$ in the range $[s_1; s_2]$ is linearly unstable, leading to a heterogeneous distribution of dry regions with humid regions. This aspect is studied in detail in the next Section.

In unsaturated porous media, the presence of the water/air interface implies the existence of a free surface energy. The pressure term $T\Delta s$ was added in [7, 15] in the