# Comptes Rendus
Mathématique

Nathan Shourick, Ibrahim Cheddadi and Pierre Saramito

A new framework for shallow approximations of incompressible flows

Volume 361 (2023), p. 1767-1783

Published online: 21 December 2023

https://doi.org/10.5802/crmath.526

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

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
In the literature, two main cases are considered: a slip-friction together with a no-penetration condition [4, 5, 7, 8, 10] and a no-slip condition [1, 9].

Our aim is to propose a general framework for facilitating the formal derivation of the thin-layer approximation of incompressible fluids, either Newtonian or non-Newtonian, with slip-friction and no-penetration boundary condition on a flat substrate. It allows both to avoid redundant mathematical hypotheses and to dramatically reduce most tedious formal calculations. While we focus on fluid mechanics applications, the present framework could be suitable to a wider class of partial differential equations from mathematical physics. The idea to provide some kind of generalization has already been suggested in the literature. In 2004, Bouchut and Westdickenberg [6] studied general varying topographies and, in 2016, Bouchut and Boyaval [5] proposed for the first time a unified framework for deriving thin-layer models for shallow free-surface flows driven by gravity, with slip-friction and no-penetration boundary condition, and under the motion by slices hypothesis. The present framework requires a minimal number of assumptions. In particular, we neither assume that the system is driven by gravity nor that the motion by slices hypothesis is satisfied.

The viscoelastic Maxwell model is studied as a proof of concept for the proposed framework. It writes [11]:

$$ \lambda \left( \frac{\partial \boldsymbol{\tau}}{\partial t} + (\boldsymbol{v} \cdot \nabla) \boldsymbol{\tau} - \nabla \boldsymbol{v} \cdot \boldsymbol{\tau} - \boldsymbol{\tau} \cdot \nabla \boldsymbol{v}^{\top} \right) + \boldsymbol{\tau} = 2\eta D(\boldsymbol{v}), \quad (1a) $$

$$ \rho \left( \frac{\partial \boldsymbol{v}}{\partial t} + (\boldsymbol{v} \cdot \nabla) \boldsymbol{v} \right) - \operatorname{div} \boldsymbol{\sigma} = \rho \boldsymbol{g}, \quad (1b) $$

$$ \operatorname{div} \boldsymbol{v} = 0, \quad (1c) $$

where the total stress tensor is $\boldsymbol{\sigma} = \boldsymbol{\tau} - p\boldsymbol{I}$, and the three unknowns are $\boldsymbol{\tau}$, the elastic stress tensor, $\boldsymbol{v}$, the velocity field and $p$, the pressure. Here, $\rho$ is the density of the fluid, assumed to be constant, and $\boldsymbol{g}$ is the constant gravity vector. The parameter $\eta$ is a viscosity and $\lambda$ is a relaxation time. The notation $D(\boldsymbol{v})$ stands for the rate of deformation tensor, that is the symmetric part of the velocity gradient, defined by $D(\boldsymbol{v}) := (\nabla \boldsymbol{v} + \nabla \boldsymbol{v}^{\top})/2$. Note that when $\lambda = 0$, the Maxwell model reduces to the Navier-Stokes one. Rather than directly performing the asymptotic analysis of the system (1), we propose in this paper to treat the following general equation:

$$ \mathcal{A}(u) - \operatorname{div} j = 0, \quad (2) $$

where $\mathcal{A}$ is a given scalar operator depending on the scalar unknown $u$ and $j$ is the equivalent of a flux. This equation is scalar but it can correspond component by component to vectorial conservation equations such as (1b), and tensorial evolution equations such as (1a), and, in general, to a large class of equations from mathematical physics.

After a presentation of the general mathematical setting in Section 2, the asymptotic analysis of (2) will be presented along with sufficient conditions for the motion by slices. This section ends with the asymptotic analysis of the general objective Gordon-Schowalter derivatives of both vector and tensor fields. Finally, Section 4 turns to applications to Navier-Stokes and Maxwell systems. To the best of our knowledge, the tools presented here are new and we hope that they will (i) help to clarify the mathematical hypothesis required for asymptotic analysis, and (ii) make formal calculations easier for future applications in the study of complex materials.

## 2. Mathematical setting

### 2.1. Geometry

We consider here a thin-layer fluid domain $\mathcal{B} \subset \mathbb{R}^3$ moving at velocity $\boldsymbol{v} : \mathcal{B} \times \mathbb{R}_+ \to \mathbb{R}^3$, with a free surface characterized by its height $h : \mathcal{B} \times \mathbb{R}_+ \to \mathbb{R}_+$, assumed small compared to the flow length.
Let $\Sigma_s$ be the three-dimensional rigid substrate on which the fluid flows. It is defined as the set of points $\mathbf{x}$ with $x_3 = 0$, and $\Omega \subset \mathbb{R}^2$ its two-dimensional counterpart such that $\Sigma_s := \Omega \times \{0\}$, see Figure 1. The free surface is defined as the set $\Sigma_f(t) := \Omega \times \{h(t)\}$. The thin-layer approximation aims at obtaining a two-dimensional problem defined over $\Omega$.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Asymptotic analysis of a free-surface flow of a shallow fluid.</figcaption></figure>

## 2.2. Evolution equations

A large part of the equations derived from mechanics can be written as one of the following equations:

$$ \lambda \frac{D\varphi}{Dt} + b(\varphi) - \operatorname{div} \mathbf{j} = f, \qquad (3a) $$

$$ \lambda \frac{\mathcal{D}_a \mathbf{u}}{\mathcal{D}t} + \mathbf{b}(\mathbf{u}) - \operatorname{div} \boldsymbol{\sigma} = \mathbf{f}, \qquad (3b) $$

$$ \text{or} \quad \lambda \frac{\mathcal{D}_a \boldsymbol{\tau}}{\mathcal{D}t} + \boldsymbol{\beta}(\boldsymbol{\tau}) - \operatorname{div} \mathbf{J} = \boldsymbol{\phi}, \qquad (3c) $$

where $\lambda$ is a time relaxation, $b$ (resp. $\mathbf{b}$ and $\boldsymbol{\beta}$) is a generic scalar-valued (resp. vector and symmetric tensor) source term, $j$ (resp. $\boldsymbol{\sigma}$ and $\mathbf{J}$) is a vector-valued (resp. tensor-valued and order 3 tensor-valued) flux and $f$ (resp. $\mathbf{f}$ and $\boldsymbol{\phi}$) is a scalar-valued (resp. vector and symmetric tensor) source term independent of the scalar (resp. vector and symmetric tensor) unknown $\varphi$ (resp. $\mathbf{u}$ and $\boldsymbol{\tau}$). Three time derivatives have also been introduced. The first one is the Lagrangian derivative $\mathcal{D}\Box/\mathcal{D}t := \partial\Box/\partial t + (\mathbf{v} \cdot \nabla)\Box$. The two others are the vector and the tensor Gordon-Schowalter derivatives, parameterized by $a \in \mathbb{R}$ and defined by

$$ \frac{\mathcal{D}_a \mathbf{u}}{\mathcal{D}t} := \frac{D\mathbf{u}}{Dt} - (W(\mathbf{v}) + aD(\mathbf{v})) \cdot \mathbf{u}, \\ \frac{\mathcal{D}_a \boldsymbol{\tau}}{\mathcal{D}t} := \frac{D\boldsymbol{\tau}}{Dt} - 2 \operatorname{sym}((W(\mathbf{v}) + aD(\mathbf{v})) \cdot \boldsymbol{\tau}), $$

where $\operatorname{sym}(\cdot)$ is the tensor symmetric part operator and $W(\mathbf{v}) := (\nabla \mathbf{v} - \nabla \mathbf{v}^\top)/2$ is the vorticity tensor, i.e. the skew-symmetric part of the velocity gradient. Note that these time derivatives are defined for a general velocity field $\mathbf{v}$, not necessarily divergence-free.

Note that all the three variants of (3) reduce component-by-component to the general case (2) for a suitable choice of the $\mathcal{A}$ that includes both the time derivative, the source and the right-hand-side terms. This approach leads to split the difficulties for the asymptotic analysis of the next Section 3: we first study (2) with a general operator $\mathcal{A}$ before we deal with the complicated terms involved by Gordon-Schowalter derivatives.
## 2.3. Boundary conditions

Only the boundary conditions satisfied on the free surface and on the substrate are involved by the asymptotic analysis. Indeed, conditions on vertical boundary walls have no effects on the reduced models.

### 2.3.1. On the free surface

For fluid flow systems such as Navier-Stokes or Maxwell (1), the boundary conditions on the free surface $\Sigma_f(t)$ write [5, 10]:

$$ \boldsymbol{\sigma} \cdot \boldsymbol{n} - \gamma \kappa \boldsymbol{n} = \mathbf{0}, \qquad (4a) $$

$$ \text{and} \quad \frac{\partial h}{\partial t} + v_1 \frac{\partial h}{\partial x_1} + v_2 \frac{\partial h}{\partial x_2} = v_3 \qquad (4b) $$

where $\boldsymbol{n}$ is the unit outward normal vector, along the direction $(-\partial h/\partial x_1, -\partial h/\partial x_2, 1)$. Eqn. (4a) expresses the normal-stress continuity condition where $\gamma$ is the surface tension at the air/fluid interface and $\kappa := -\text{div} \boldsymbol{n}$ is the local mean curvature while (4b) is the usual free surface kinematic equation. Note that (4a) is a non-homogeneous Neumann boundary condition. By extension, for our general scalar equation (2), we consider the general condition on free surface $\Sigma_f(t)$:

$$ \boldsymbol{j} \cdot \boldsymbol{n} = q. \qquad (5) $$

### 2.3.2. On the substrate

On the substrate $\Sigma_s$, we consider:

$$ \boldsymbol{v} \cdot \boldsymbol{n} = 0, \qquad (6a) $$

$$ \boldsymbol{\sigma}_{nt} = k \boldsymbol{v}_t, \qquad (6b) $$

where $\boldsymbol{\sigma}_{nt} := \boldsymbol{\sigma} \cdot \boldsymbol{n} - (\boldsymbol{n} \cdot \boldsymbol{\sigma} \cdot \boldsymbol{n})\boldsymbol{n}$ is the tangential part of the normal stress, $\boldsymbol{v}_t := \boldsymbol{v} - (\boldsymbol{v} \cdot \boldsymbol{n})\boldsymbol{n}$ is the tangential part of the velocity. Relation (6a) expresses the *no-penetration* of the flow material across the substrate surface while (6b) is the tangential *slip-friction* condition, where $k$ is the friction coefficient, which is constant for laminar flows.

Condition (6) is the usual starting point of most *motion by slices* approximate models. Note that a popular alternative is the *no-slip* condition $\boldsymbol{v} = \mathbf{0}$, which is the starting point of the so-called *lubrication* approximate models [1, 2, 9]. For simplicity, the later will not be considered and we assume (6) all along this paper.

# 3. The framework

## 3.1. Notations

### 3.1.1. Asymptotic analysis

Let $H$, $L$ and $U$ be the characteristic height, length and velocity, respectively. For instance $L$ is the substrate length. The asymptotic analysis consists in exploiting the small, yet fixed, aspect ratio

$$ \epsilon := H/L $$

in order to reduce the three-dimensional problem to a bidimensional one. For this purpose, we are going to scale the spatial coordinates with this parameter, distinguishing the contribution from planar and horizontal dimensions. We then assume each *dimensionless* unknown admits a Taylor expansion with respect to $\epsilon$ at second order, which, to fix ideas, writes

$$ \varphi = \varphi^{(0)} + \varphi^{(1)}\epsilon + \varphi^{(2)}\epsilon^2 + O(\epsilon^3). \qquad (7) $$
By convention, we consider that both the height and the $x_3$-coordinate are exactly of order one while both the $x_1$-coordinate and $x_2$-coordinate do not depend on $\epsilon$:

$$ \mathbf{x}_s = L\mathbf{x}_s^{(0)}, \quad x_3 = Lx_3^{(1)}\epsilon, \quad t = \frac{L}{U}t^{(0)}, \qquad (8) $$

where the s index refers to the planar components of a given variable, e.g. $\mathbf{x}_s := (x_1, x_2)$. We deduce the following scaling for differential operators:

$$ \nabla_s = \frac{1}{L} \nabla_s^{(0)}, \quad \frac{\partial}{\partial x_3} = \frac{1}{\epsilon L} \frac{\partial}{\partial x_3^{(1)}}, \quad \frac{\partial}{\partial t} = \frac{U}{L} \frac{\partial}{\partial t^{(0)}}, \qquad (9) $$

where $\nabla_s := (\partial/\partial x_1, \partial/\partial x_2)^\top$ is the planar gradient. We also define the planar divergence operator $\text{div}_s \square := \partial\square_{i_1...i_{p-1}}/\partial x_1 + \partial\square_{i_1...i_{p-12}}/\partial x_2$, where $p$ is the order of the given tensor, and the planar Laplacian $\Delta_s := \nabla_s \cdot \nabla_s$. Similarly, we assume that the dimensionless velocity field $v/U$ and the dimensionless height field $h/L$ satisfy an expansion of the form (7) with $v_3^{(0)} = 0$ and $h^{(0)} = 0$. This choice is justified by our wish to have $\epsilon U$ as the characteristic vertical speed and $H = \epsilon L$ as the characteristic height.

Observe that the notation $\varphi^{(k)}$ makes it possible to immediately recognise the order with respect to $\epsilon$ of the considered term. For instance, the product $\varphi^{(k)}\psi^{(l)}$ will be associated to a term of order $O(\epsilon^{k+l})$, while the quotient $\varphi^{(k)}/\psi^{(l)}$ will be associated to a term of order $O(\epsilon^{k-l})$.

From now until Section 4, we assume that every field, parameter and differential operator is dimensionless, for convenience. In particular, when we refer to the dimensionless height (resp. vertical velocity), we mean the height scaled by $L$ (resp. $U$), in order to keep the property presented in the previous paragraph.

### 3.1.2. Depth averaging

Depth averaging is the final stage of the reduction of equations. It is defined through the following operator:

**Definition 1 (depth average).** Let $\varphi : \mathcal{B} \times \mathbb{R}_+ \to \mathbb{R}$ a continuous scalar field with respect to $x_3$. Its depth average is defined by the quantity

$$ \langle \varphi \rangle_h(\mathbf{x}_s, t) := \begin{cases} \frac{1}{h} \int_0^h \varphi((\mathbf{x}_s, x_3), t) dx_3, & \forall (\mathbf{x}_s, t) \in \Omega \times \mathbb{R}_+ \quad \text{when } h(\mathbf{x}_s, t) \neq 0, \\ \varphi((\mathbf{x}_s, 0), t) & \text{otherwise.} \end{cases} $$

This linear operator extends componentwise to any vector or tensor field. We also define the depth variance and covariance:

**Definition 2 (depth variance and covariance).** Let $\varphi$ and $\psi$ be two scalar fields defined in $\mathcal{B} \times \mathbb{R}_+$, continuous with respect to $x_3$.

* The depth variance of $\varphi$ is defined by the quantity

$$ \text{var}_h(\varphi) := \langle (\varphi - \langle \varphi \rangle_h)^2 \rangle_h. $$

* The depth covariance of $\varphi$ and $\psi$ is defined by the quantity

$$ \text{cov}_h(\varphi, \psi) := \langle (\varphi - \langle \varphi \rangle_h)(\psi - \langle \psi \rangle_h) \rangle_h. $$

In particular, $\text{cov}_h(\varphi, \varphi) = \text{var}_h(\varphi)$.

Those operators satisfy the Cauchy–Schwartz inequality

$$ |\text{cov}_h(\varphi, \psi)|^2 \le \text{var}_h(\varphi)\text{var}_h(\psi), \qquad (10a) $$
and the König-Huygens theorem-like

$$ \mathrm{var}_h(\varphi) = \langle\varphi^2\rangle_h - \langle\varphi\rangle_h^2, \qquad (10b) $$

$$ \mathrm{cov}_h(\varphi, \psi) = \langle\varphi\psi\rangle_h - \langle\varphi\rangle_h \langle\psi\rangle_h. \qquad (10c) $$

Note the following useful identities:

$$ \langle x_3 \rangle_h = h/2 \quad \text{and} \quad \mathrm{var}_h x_3 = h^2/12. \qquad (11) $$

## 3.2. Analysis of the general scalar problem

Considering equation (2) together with the Neumann condition (5), let us perform the asymptotic analysis of the following general scalar problem:

**Problem 3 (general scalar).** Find $u$, defined in $\mathcal{B} \times \mathbb{R}_+$, such that for any time $t \in \mathbb{R}_+$

$$ \mathcal{A}(u(t)) - \mathrm{div} \mathbf{j}(t) = \mathbf{0} \quad \text{in } \mathcal{B}(t), \qquad (12a) $$

$$ \mathbf{j}(t) \cdot \mathbf{n} = q \quad \text{on } \Sigma_f(t). \qquad (12b) $$

Recall that this scalar problem extends componentwise to a vector or a tensor one. Assume that $u$ is a solution of (12) and that the following expansions hold:

$$ \mathcal{A} = \mathcal{A}^{(-1)} \frac{1}{\epsilon} + \mathcal{A}^{(0)} + \mathcal{A}^{(1)} \epsilon + O(\epsilon^2), \qquad (13a) $$

$$ \mathbf{j} = \mathbf{j}^{(-1)} \frac{1}{\epsilon} + \mathbf{j}^{(0)} + \mathbf{j}^{(1)} \epsilon + \mathbf{j}^{(2)} \epsilon^2 + O(\epsilon^3), \qquad (13b) $$

$$ q = q^{(-1)} \frac{1}{\epsilon} + q^{(0)} + q^{(1)} \epsilon + q^{(2)} \epsilon^2 + O(\epsilon^3) \qquad (13c) $$

Then (12a) leads to

$$ -\frac{\partial j_3^{(-1)}}{\partial x_3^{(1)}} \frac{1}{\epsilon^2} + \sum_{k=-1}^{1} \left( \mathcal{A}^{(k)} - \mathrm{div}_s \mathbf{j}_s^{(k)} - \frac{\partial j_3^{(k+1)}}{\partial x_3^{(1)}} \right) \epsilon^k + O(\epsilon^2) = 0, $$

while the boundary condition (12b) becomes

$$ (j_3^{(-1)} - q^{(-1)}) \frac{1}{\epsilon} + \sum_{k=-1}^{1} (j_3^{(k+1)} - j_s^{(k)} \cdot \nabla_s h^{(1)} - q^{(k+1)}) \epsilon^{k+1} + O(\epsilon^2) = 0 \quad \text{on } \Sigma_f(t). $$

By identification, we get, for any $k \in \{-1, 0, 1\}$

$$ \frac{\partial j_3^{(-1)}}{\partial x_3^{(1)}} = 0, \qquad (14a) $$

$$ \mathcal{A}^{(k)} - \mathrm{div}_s \mathbf{j}_s^{(k)} = \frac{\partial j_3^{(k+1)}}{\partial x_3^{(1)}}, \qquad (14b) $$

$$ j_3^{(-1)} = q^{(-1)} \qquad \text{on } \Sigma_f(t), \qquad (14c) $$

$$ j_3^{(k+1)} = \nabla_s h^{(1)} \cdot \mathbf{j}_s^{(k)} + q^{(k+1)} \quad \text{on } \Sigma_f(t). \qquad (14d) $$

From (14a) and (14c), we also deduce $j_3^{(-1)} = q^{(-1)}$ in $\Omega \times \mathbb{R}_+$.

**Definition 4 (hypotheses $\mathcal{H}^{(k)}$).** For any $k \in \{-1, 0, 1\}$, we say hypothesis $\mathcal{H}^{(k)}$ holds if and only if $\mathcal{A}^{(k)} = 0$ and $\mathbf{j}_s^{(k)} = \mathbf{0}$.
**Lemma 5 (vertical flux and Neumann BC).** Assume expansions (13) are satisfied. Then,

* the term of order $-1$ of the vertical flux is given by

$$j_3^{(-1)} = q^{(-1)} \quad \text{in } \Omega \times \mathbb{R}_+ \qquad (15a)$$

* if hypothesis $\mathcal{H}^{(k)}$ also holds, for $k \in \{-1, 0, 1\}$,

$$j_3^{(k+1)} = q^{(k+1)} \quad \text{in } \Omega \times \mathbb{R}_+ \qquad (15b)$$

**Proof.** The first point has already been proven in the main text. Under hypothesis $\mathcal{H}^{(k)}$, (14b) reduces to

$$\frac{\partial j_3^{(k+1)}}{\partial x_3^{(1)}} = 0 \quad \text{in } \Omega \times \mathbb{R}_+,$$

whence the result, from the boundary condition (14d). $\square$

This lemma is the starting point for obtaining the property of motion by slices (see Corollary 8).

It remains to average in depth the problem defined by equations (14a)–(14b) together with boundary conditions (14c)–(14d). This can be performed thanks to Corollary 26 in the appendix, which gives the

**Problem 6 (reduced general scalar).** Find $\langle u^{(0)} \rangle_h$ and $\langle u^{(1)} \rangle_h$, defined in $\Omega \times \mathbb{R}_+$, such that for any time $t \in \mathbb{R}_+$ and any index $k \in \{-1, 0, 1\}$

$$h^{(1)} \langle \mathcal{A}^{(k)} \rangle_h - \operatorname{div}_s \left( h^{(1)} \langle j_s^{(k)} \rangle_h \right) = q^{(k+1)} - j_3^{(k+1)} (x_3 = 0). \qquad (16)$$

**Theorem 7 (asymptotic analysis of the general problem).** Let $u$ be the solution of Problem 3, which is assumed to be satisfying expansion (13), and $\langle u^{(0)} \rangle_h$ and $\langle u^{(1)} \rangle_h$ be solutions of Problem 6. Then,

$$\langle u \rangle_h = \langle u^{(0)} \rangle_h + \langle u^{(1)} \rangle_h \epsilon + O(\epsilon^2). \qquad (17)$$

In other words, the reconstructed solution $\langle u^{(0)} \rangle_h + \langle u^{(1)} \rangle_h \epsilon$ is an approximation of order $O(\epsilon^2)$ of the averaged exact solution $\langle u \rangle_h$.

**Corollary 8 (motion by slices).** Under the same assumptions than the above Theorem 7, let us assume in addition that $\mathbf{j} = v\nabla u + \mathbf{m}$, where $v > 0$ and $\mathbf{m}$ is a vector field that does not depend on $\nabla u$. If $\mathbf{m}$ expands as $\mathbf{m} = \mathbf{m}^{(0)} + \mathbf{m}^{(1)}\epsilon + \mathbf{O}(\epsilon^2)$, $q^{(-1)} = 0$ and $v = O(1)$, then

$$\frac{\partial u^{(0)}}{\partial x_3^{(1)}} = 0.$$

**Proof.** By construction, $j_3^{(-1)} = v\partial u^{(0)}/\partial x_3^{(1)}$. The result follows then from (15a). $\square$

**Remark 9 (no-slip condition).** More generally, if $\mathbf{m}$ expands as $\mathbf{m} = \mathbf{m}^{(-1)} / \epsilon + \mathbf{m}^{(0)} + \mathbf{m}^{(1)}\epsilon + \mathbf{m}^{(2)}\epsilon^2 + \mathbf{O}(\epsilon^3)$ and $v = O(\epsilon^n)$, with $n \in \{0, \dots, 3\}$, then

$$\frac{\partial u^{(0)}}{\partial x_3^{(1)}} = \frac{1}{v} (q^{(n-1)} - m_3^{(n-1)}), \quad \text{that is} \quad \langle u^{(0)} \rangle_h = u^{(0)}(x_3 = 0) + \frac{h^{(1)}}{2v} q^{(n-1)} - \left\langle \int_0^{x_3} \frac{1}{v} m_3^{(n-1)} dx_3 \right\rangle_h. \qquad (18)$$

Instead of the Neumann boundary condition (12b), we could consider a Dirichlet condition on $u$. In that case, we do no more have information on the normal flux $\mathbf{j} \cdot \mathbf{n}$ on the substrate. This could be a difficulty since this term is involved in the right-hand-side of (16). Note that a popular case of such Dirichlet condition is the no-slip condition on the substrate, which the starting point of lubrication shallow models. In that case, there is no more motion by slices, but this asymptotic expansion still allows to obtain an explicit expression of the zeroth order depth average of $u$.
Moreover, from (15b) and under assumption $\mathcal{H}^{(k)}$ (see Lemma 5), it is possible to extend the result to any $k \in \{-1, 0, 1\}$ as

$$ \langle u^{(k+1)} \rangle_h = u^{(k+1)}(x_3 = 0) + \frac{h^{(1)}}{2v} q^{(n+k)} - \left\langle \int_0^{x_3} \frac{1}{v} m_3^{(n+k)} dx_3 \right\rangle_h, \quad (19) $$

See e.g. [1, §3.1] or [2] for practical examples of lubrication problems.

**Corollary 10 (asymptotic analysis of the pressure).** For any $i \in \{1, 2, 3\}$, assume the $i$-th component of a vector field $\mathbf{u}$ is solution of a problem of the form (12) with vector flux $(\boldsymbol{\sigma}_{ij})_{1 \le j \le 3}$ and operator $\mathcal{A}_i$. Assume in addition that these three problems satisfy the expansion (13). Assume also that the tensor flux writes $\boldsymbol{\sigma} = \boldsymbol{\tau} - p\mathbf{I}$, where $\boldsymbol{\tau}$ is a tensor field and $p$ is a scalar field, called the pressure, both of them being of order at most $O(1/e)$. Then the depth average of the pressure $p$ is given by

$$ \begin{aligned} \langle p^{(-1)} \rangle_h &= \langle \tau_{33}^{(-1)} \rangle_h - q^{(-1)}, \\ \langle p^{(k+1)} \rangle_h &= \langle \tau_{33}^{(k+1)} \rangle_h - q^{(k+1)} - \left\langle \int_{x_3}^h \left( \operatorname{div}_s \boldsymbol{\tau}_{3s}^{(k)} - \mathcal{A}_3^{(k)} \right) dx_3 \right\rangle_h - \boldsymbol{\tau}_{3s}^{(k)}(x_3 = h) \cdot \nabla_s h^{(1)}, \end{aligned} $$ 
for any $k \in \{-1, 0, 1\}$.

**Proof.** The first result is immediate from (15a). From (14b), we get

$$ \frac{\partial p^{(k+1)}}{\partial x_3} = \frac{\partial \tau_{33}^{(k+1)}}{\partial x_3} + \operatorname{div}_s \boldsymbol{\tau}_{3s}^{(k)} - \mathcal{A}_3^{(k)}, $$

and then, by integrating between $x_3$ and $h$

$$ p^{(k+1)}(x_3 = h) - p^{(k+1)} = \tau_{33}^{(k+1)}(x_3 = h) - \tau_{33}^{(k+1)} + \int_{x_3}^{h} \left( \operatorname{div}_s \boldsymbol{\tau}_{3s}^{(k)} - \mathcal{A}_3^{(k)} \right) dx_3. $$

From the Neumann boundary condition (14d), we obtain the relation

$$ \tau_{33}^{(k+1)}(x_3 = h) - p^{(k+1)}(x_3 = h) = \boldsymbol{\tau}_{3s}^{(k)}(x_3 = h) \cdot \nabla_s h^{(1)} + q^{(k+1)}. $$

A depth integration allows then to conclude. $\square$

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

**Remark 11 (simplifications).** If $\boldsymbol{\tau}_{3s}^{(k)}$ is independent upon $x_3$, we have the much simpler result

$$ \langle p^{(k+1)} \rangle_h = \langle \tau_{33}^{(k+1)} \rangle_h - q^{(k+1)} - \operatorname{div}_s \left( h^{(1)} \boldsymbol{\tau}_{3s}^{(k)} \right) + \frac{h^{(1)}}{2} \operatorname{div}_s \boldsymbol{\tau}_{3s}^{(k)} + \left\langle \int_{x_3}^{h} \mathcal{A}_3^{(k)} dx_3 \right\rangle_h. $$

Similarly, if $\mathcal{A}_3^{(k)}$ is also independent of $x_3$, we have the even simpler result

$$ \langle p^{(k+1)} \rangle_h = \langle \tau_{33}^{(k+1)} \rangle_h - q^{(k+1)} - \operatorname{div}_s \left( h^{(1)} \boldsymbol{\tau}_{3s}^{(k)} \right) + \frac{h^{(1)}}{2} \operatorname{div}_s \boldsymbol{\tau}_{3s}^{(k)} + \frac{h^{(1)}}{2} \mathcal{A}_3^{(k)}. $$

### 3.3. Time derivatives

The aim of this section is to derive the depth average of the time derivatives that we introduced in Section 2.2, namely the material derivative $D/Dt$ and the vector and tensor Gordon-Schowalter derivatives $\mathcal{D}_a/\mathcal{D}t$, and express them as functions of their equivalents in the plane $\Omega$. We will do so under the assumptions that lead to motion by slices (see Corollary 8):

$$ \frac{\partial \boldsymbol{v}_s^{(0)}}{\partial x_3^{(1)}} = \mathbf{0}, \qquad (20) $$

and the fact that the velocity field $\boldsymbol{v}$ is divergence-free (incompressibility) and satisfies to the no-penetration constraint on the substrate. To exploit this, we will first need to perform the asymptotic analysis of the divergence-free condition.
### 3.3.1. Divergence-free condition

Let us start by the asymptotic analysis of the divergence-free condition

$$ \operatorname{div} \boldsymbol{v} = 0, \qquad (21) $$

under the boundary conditions

$$ \boldsymbol{v} \cdot \boldsymbol{n} = 0 \iff v_3 = 0 \quad \text{on } \Sigma_s, \qquad (22a) $$

$$ \frac{D_s h}{Dt} = v_3 \quad \text{on } \Sigma_f(t), \qquad (22b) $$

where we noted $\frac{D_s}{Dt} := \frac{\partial}{\partial t} + (\boldsymbol{v}_s \cdot \nabla_s) \boldsymbol{v}_s$. The material derivative advected by the velocity $\boldsymbol{v}_s$.

**Proposition 12 (averaged vertical velocity).** Assume that $\boldsymbol{v}$ is divergence-free and that the property of motion by slices (20) and the boundary conditions (22) are satisfied. Then the vertical velocity $v_3$ is given by:

$$ v_3 = -\operatorname{div}_s (\boldsymbol{v}_s^{(0)}) x_3 - \langle \operatorname{div}_s (\boldsymbol{v}_s^{(1)}) \rangle_{x_3} \epsilon + O(\epsilon^3), \quad \text{that is } v_3^{(1)} = -\operatorname{div}_s (\boldsymbol{v}_s^{(0)}) x_3^{(1)}. \qquad (23) $$

**Proof.** Note that $\operatorname{div} \boldsymbol{v} = \operatorname{div}_s \boldsymbol{v}_s + \partial v_3 / \partial x_3 = 0$. Therefore, the result follows directly from an integration of the equation (21) using the hypothesis of motion by slices and the boundary condition (22a). $\square$

**Corollary 13 (averaged free surface evolution).** Under the same conditions as in the above proposition, the height $h$ is solution of the equation

$$ \frac{\partial h}{\partial t} + \operatorname{div}_s (h \boldsymbol{v}_s^{(0)}) + \operatorname{div}_s (h \langle \boldsymbol{v}_s^{(1)} \rangle_h) \epsilon = O(\epsilon^3). \qquad (24) $$

**Proof.** Immediate from Proposition 12, Corollary 26 (see Appendix A) and boundary condition (22b). $\square$

By identification, we immediately obtain the equation satisfied by the height at order $O(\epsilon)$:

$$ \frac{\partial h^{(1)}}{\partial t} + \operatorname{div}_s (h^{(1)} \boldsymbol{v}_s^{(0)}) = 0. \qquad (25) $$

Solving (25) instead of (24) leads to an error of order $O(\epsilon)$ since<fnref n="1" /> $h/\epsilon = h^{(1)} + O(\epsilon)$.

### 3.3.2. Lagrangian derivative

The depth average of the Lagrangian derivative is expressed as a function of the planar material derivative, advected by the velocity $\boldsymbol{v}_s^{(0)}$, noted

$$ \frac{D_s^{(0)} \varphi}{Dt} := \frac{\partial \varphi}{\partial t} + (\boldsymbol{v}_s^{(0)} \cdot \nabla_s) \varphi, $$

for any scalar field $\varphi$ defined in $\Omega \times \mathbb{R}_+$. It extends componentwise to any vector or tensor field. Thanks to the new mass conservation law (24), we deduce the

**Proposition 14 (averaged Lagrangian derivative).** Assume that $\boldsymbol{v}$ is divergence-free and that the property of motion by slices (20) and the boundary conditions (22) are satisfied. Then, any scalar field $u$ defined in $\mathcal{B} \times \mathbb{R}_+$ satisfies the relation

$$ \left\langle \frac{Du}{Dt} \right\rangle_h = \frac{D_s^{(0)} \langle u \rangle_h}{Dt} + O(\epsilon), \qquad (26) $$

<footnote n="1">¹As a reminder, the height was scaled by the characteristic length $L$, not by the characteristic height $H$, as mentioned in Section 3.1.1.</footnote>
Moreover, if we assume $\operatorname{cov}_h(\nabla_s u, \mathbf{v}_s^{(1)}) = O(\epsilon)$, then we have a better approximation:

$$ \left\langle \frac{Du}{Dt} \right\rangle_h = \frac{D_s^{(0)} \langle u \rangle_h}{Dt} + (\langle \mathbf{v}_s^{(1)} \rangle_h \cdot \nabla_s) \langle u \rangle_h \epsilon + O(\epsilon^2), \quad (27) $$

**Proof.** We show this property in the second case, the first result being independent of the assumption $\operatorname{cov}_h(\nabla_s u, \mathbf{v}_s^{(1)}) = O(\epsilon)$. The calculation is performed in two steps, by noting that $D/Dt = D_s/Dt + v_3 \partial/\partial x_3 = D_s/Dt + v_3^{(1)} \partial/\partial x_3^{(1)} + O(\epsilon)$ and by using the linearity of the operator $\langle \cdot \rangle_h$.

(1) Two direct applications of the Leibniz formula (Corollary 25 in Appendix A) and one of the formula (10c) gives

$$ \begin{aligned} h \left\langle \frac{D_s u}{Dt} \right\rangle_h &= \frac{D_s^{(0)}}{Dt} (h \langle u \rangle_h) - u(x_3 = h) \frac{D_s^{(0)} h}{Dt} + h \langle (\mathbf{v}_s^{(1)} \cdot \nabla_s) u \rangle_h \epsilon + O(\epsilon^3) \\ &= h \left( \frac{D_s^{(0)} \langle u \rangle_h}{Dt} + (\langle \mathbf{v}_s^{(1)} \rangle_h \cdot \nabla_s) \langle u \rangle_h \epsilon \right) + h \operatorname{cov}_h(\nabla_s u, \mathbf{v}_s^{(1)}) \epsilon \\ &\quad + (\langle u \rangle_h - u(x_3 = h)) \left( \frac{D_s^{(0)} h}{Dt} + (\langle \mathbf{v}_s^{(1)} \rangle_h \cdot \nabla_s) h \epsilon \right) + O(\epsilon^3). \end{aligned} $$

(2) A direct application of the appendix formula (47e), after expansion of the vertical velocity using the Proposition 12 and use of the hypothesis of motion by slices, gives

$$ \begin{aligned} h \left\langle v_3 \frac{\partial u}{\partial x_3} \right\rangle_h &= - (h \operatorname{div}_s \mathbf{v}_s^{(0)} + h \langle \operatorname{div}_s \mathbf{v}_s^{(1)} \rangle_h \epsilon) \left\langle x_3 \frac{\partial u}{\partial x_3} \right\rangle_h + O(\epsilon^3) \\ &= (h \operatorname{div}_s \mathbf{v}_s^{(0)} + h \langle \operatorname{div}_s \mathbf{v}_s^{(1)} \rangle_h \epsilon) (\langle u \rangle_h - u(x_3 = h)) + O(\epsilon^3). \end{aligned} $$

The sum of the two above results and the application of equality (24) allow to conclude. $\square$

This result still holds true when replacing $u$ with a vector or a tensor field, in particular with the velocity field $\mathbf{v}$.

### 3.3.3. Gordon–Schowalter derivative of a vector

First, the depth average of the Gordon–Schowalter derivative of a vector is expressed as a function of its plane equivalent, with respect to the zeroth-order planar velocity $\mathbf{v}_s^{(0)}$:

$$ \frac{\mathcal{D}_a^{(0)} \mathbf{w}_s}{\mathcal{D}t} := \frac{D_s^{(0)} \mathbf{w}_s}{Dt} - W_s(\mathbf{v}_s^{(0)}) \cdot \mathbf{w}_s - aD_s(\mathbf{v}_s^{(0)}) \cdot \mathbf{w}_s, \quad (28) $$

for any planar vector $\mathbf{w}_s$ defined in $\Omega \times \mathbb{R}_+$.

**Theorem 15 (Gordon–Schowalter derivative of a vector).** Let $\mathbf{u}$ be a vector field defined in $\mathcal{B} \times \mathbb{R}_+$ such that $\operatorname{var}_h(\mathbf{u}) = O(\epsilon)$ (closure assumption). Assume that $\mathbf{v}$ is divergence-free and that the property of motion by slices (20) and the boundary conditions (22) are satisfied. Then the depth-average of the Gordon–Schowalter derivative of $\mathbf{u}$ is given by

$$ \left( \left\langle \frac{\mathcal{D}_a \mathbf{u}}{\mathcal{D}t} \right\rangle_h \right)_s = \frac{\mathcal{D}_a^{(0)} \langle \mathbf{u}_s \rangle_h}{\mathcal{D}t} - \frac{1+a}{2} \left\langle \frac{\partial \mathbf{v}_s^{(1)}}{\partial x_3^{(1)}} u_3 \right\rangle_h - \left( \langle [W_s(\mathbf{v}_s^{(1)}) + aD_s(\mathbf{v}_s^{(1)})] \cdot \mathbf{u}_s \rangle_h \right. \\ \left. + \frac{1+a}{2} \left\langle \frac{\partial \mathbf{v}_s^{(2)}}{\partial x_3^{(1)}} u_3 \right\rangle_h + \frac{(1-a)h^{(1)}}{4} \langle u_3 \rangle_h \nabla_s \operatorname{div}_s \mathbf{v}_s^{(0)} \right) \epsilon + \mathbf{O}(\epsilon^2), \quad (29a) $$
$$ \left( \left\langle \frac{\mathcal{D}_a \mathbf{u}}{\mathcal{D}t} \right\rangle_h \right)_3 = \frac{D_s^{(0)} \langle u_3 \rangle_h}{Dt} + \frac{1-a}{2} \left\langle \frac{\partial \mathbf{v}_s^{(1)}}{\partial x_3^{(1)}} \cdot \mathbf{u}_s \right\rangle_h + a \operatorname{div}_s (\mathbf{v}_s^{(0)}) \langle u_3 \rangle_h \\ + \left( \frac{1-a}{2} \left\langle \frac{\partial \mathbf{v}_s^{(2)}}{\partial x_3^{(1)}} \cdot \mathbf{u}_s \right\rangle_h + \frac{(1+a)h^{(1)}}{4} \nabla_s \operatorname{div}_s \mathbf{v}_s^{(0)} \cdot \langle \mathbf{u}_s \rangle_h + a \langle \operatorname{div}_s (\mathbf{v}_s^{(1)}) u_3 \rangle_h \right) \epsilon + O(\epsilon^2). \quad (29b) $$

**Proof.** Recall that

$$ \frac{\mathcal{D}_a \mathbf{u}}{\mathcal{D}t} := \frac{D\mathbf{u}}{Dt} - W(\mathbf{v}) \cdot \mathbf{u} - aD(\mathbf{v}) \cdot \mathbf{u}. $$

As in the previous demonstration, we will treat each of the terms separately, and we use the previous result for the first term. The expansions (7), (8) and (9) applied to the velocity field, and the motion by slices assumption lead to

$$ 2W(\mathbf{v}) = \begin{pmatrix} 2W_s(\mathbf{v}_s^{(0)}) + 2W_s(\mathbf{v}_s^{(1)})\epsilon & \frac{\partial \mathbf{v}_s^{(1)}}{\partial x_3^{(1)}} + \left(\frac{\partial \mathbf{v}_s^{(2)}}{\partial x_3^{(1)}} - \nabla_s v_3^{(1)}\right)\epsilon \\ \left(\nabla_s v_3^{(1)} - \frac{\partial \mathbf{v}_s^{(2)}}{\partial x_3^{(1)}}\right)\epsilon - \frac{\partial \mathbf{v}_s^{(1)}}{\partial x_3^{(1)}} & 0 \end{pmatrix} + \mathbf{O}(\epsilon^2), $$

$$ 2D(\mathbf{v}) = \begin{pmatrix} 2D_s(\mathbf{v}_s^{(0)}) + 2D_s(\mathbf{v}_s^{(1)})\epsilon & \frac{\partial \mathbf{v}_s^{(1)}}{\partial x_3^{(1)}} + \left(\frac{\partial \mathbf{v}_s^{(2)}}{\partial x_3^{(1)}} + \nabla_s v_3^{(1)}\right)\epsilon \\ \left(\nabla_s v_3^{(1)} + \frac{\partial \mathbf{v}_s^{(2)}}{\partial x_3^{(1)}}\right)\epsilon + \frac{\partial \mathbf{v}_s^{(1)}}{\partial x_3^{(1)}} & 2\frac{\partial v_3^{(1)}}{\partial x_3^{(1)}} + 2\frac{\partial v_3^{(2)}}{\partial x_3^{(1)}}\epsilon \end{pmatrix} + \mathbf{O}(\epsilon^2). $$

Therefore, the second and third terms in the definition of the Gordon-Schowalter derivative are

$$ W(\mathbf{v}) \cdot \mathbf{u} = \left( \begin{aligned} & W_s(\mathbf{v}_s^{(0)}) \cdot \mathbf{u}_s + \frac{1}{2} \frac{\partial \mathbf{v}_s^{(1)}}{\partial x_3^{(1)}} u_3 + \left( W_s(\mathbf{v}_s^{(1)}) \cdot \mathbf{u}_s + \frac{\partial \mathbf{v}_s^{(2)}}{\partial x_3^{(1)}} \frac{u_3}{2} - \nabla_s v_3^{(1)} \frac{u_3}{2} \right) \epsilon + \mathbf{O}(\epsilon^2) \\ & - \frac{1}{2} \frac{\partial \mathbf{v}_s^{(1)}}{\partial x_3^{(1)}} \cdot \mathbf{u}_s + \frac{\epsilon}{2} \left( \nabla_s v_3^{(1)} - \frac{\partial \mathbf{v}_s^{(2)}}{\partial x_3^{(1)}} \right) \cdot \mathbf{u}_s + \mathbf{O}(\epsilon^2) \end{aligned} \right), $$

$$ D(\mathbf{v}) \cdot \mathbf{u} = \left( \begin{aligned} & D_s(\mathbf{v}_s^{(0)}) \cdot \mathbf{u}_s + \frac{1}{2} \frac{\partial \mathbf{v}_s^{(1)}}{\partial x_3^{(1)}} u_3 + \left( D_s(\mathbf{v}_s^{(1)}) \cdot \mathbf{u}_s + \frac{\partial \mathbf{v}_s^{(2)}}{\partial x_3^{(1)}} \frac{u_3}{2} + \nabla_s v_3^{(1)} \frac{u_3}{2} \right) \epsilon + \mathbf{O}(\epsilon^2) \\ & \frac{1}{2} \frac{\partial \mathbf{v}_s^{(1)}}{\partial x_3^{(1)}} \cdot \mathbf{u}_s + \frac{\partial \mathbf{v}_s^{(1)}}{\partial x_3^{(1)}} u_3 + \left( \frac{1}{2} \left[ \nabla_s v_3^{(1)} + \frac{\partial \mathbf{v}_s^{(2)}}{\partial x_3^{(1)}} \right] \cdot \mathbf{u}_s + \frac{\partial v_3^{(2)}}{\partial x_3^{(1)}} u_3 \right) \epsilon + \mathbf{O}(\epsilon^2) \end{aligned} \right), $$

where we used the Proposition 12. Depth averaging is easily done using the motion by slices assumption. However, it is necessary to treat the terms consisting in a product of two terms, a priori dependent on the coordinate $x_3$, carefully. We deal with the term $u_3 \nabla_s v_3^{(1)} = -x_3^{(1)} \nabla_s \operatorname{div}_s (\mathbf{v}_s^{(0)}) u_3$ for the example, whose average can be calculated through $\langle x_3^{(1)} u_3 \rangle_h$. By successively applying the formulas (10c) and (11), we get

$$ \langle x_3^{(1)} u_3 \rangle_h = \frac{h^{(1)}}{2} \langle u_3 \rangle_h + \operatorname{cov}_h(x_3^{(1)}, u_3) \approx \frac{h^{(1)}}{2} \langle u_3 \rangle_h. $$

The closure assumption allows then to prove the approximation since, from the result (10a) and the formula (11),

$$ |\operatorname{cov}_h(x_3^{(1)}, u_3)|^2 \le \operatorname{var}_h(x_3^{(1)}) \operatorname{var}_h(u_3) = \frac{h^{(1)2}}{12} \operatorname{var}_h(u_3) = O(\epsilon). $$

The linearity of $\langle \cdot \rangle_h$ and the Proposition 14 on the asymptotic analysis of the Lagrangian derivative finally lead to the conclusion. $\square$

**Remark 16 (average of a product).** The expression of the depth-average of the Gordon-Schowalter derivative involves averages of products, which is not desirable for closing the reduced problem: we would prefer instead a product of averaged quantities. The simplification of
those products is actually case-dependent: one at least of the two quantities involved by the product should be independent upon $x_3$. For instance, either the first-order planar velocity $v_s^{(1)}$ or the vector field $u$ involved by the Gordon–Schowalter derivative should be independent upon $x_3$.

### 3.3.4. Gordon–Schowalter derivative of a tensor

Then, the depth average of the Gordon–Schowalter derivative of a tensor is expressed as a function of its plane equivalent, with respect to the zeroth-order planar velocity $v_s^{(0)}$:

$$ \frac{\mathcal{D}_a^{(0)} \boldsymbol{\sigma}_{ss}}{\mathcal{D}t} := \frac{D_s^{(0)} \boldsymbol{\sigma}_{ss}}{Dt} - 2 \operatorname{sym}[W_s(\boldsymbol{v}_s^{(0)}) \cdot \boldsymbol{\sigma}_{ss} - aD_s(\boldsymbol{v}_s^{(0)}) \cdot \boldsymbol{\sigma}_{ss}], \quad (30) $$

for any planar tensor $\boldsymbol{\sigma}_{ss}$ in $\Omega \times \mathbb{R}_+$.

**Theorem 17 (Gordon–Schowalter derivative of a tensor).** Let $\boldsymbol{\tau}$ be a symmetric tensor field defined in $\mathcal{B} \times \mathbb{R}_+$, whose zeroth order vertical variation with respect to $\epsilon$ is negligible, that is $\operatorname{var}_h(\boldsymbol{\tau}^{(0)}) = O(\epsilon)$ (closure assumption). Assume that $\boldsymbol{v}$ is divergence-free and that the property of motion by slices (20) and the boundary conditions (22) are satisfied. Then the depth-average of the Gordon–Schowalter derivative of $\boldsymbol{\tau}$ is given by

$$ \left\langle \left\langle \frac{\mathcal{D}_a \boldsymbol{\tau}}{\mathcal{D}t} \right\rangle_h \right\rangle_{ss} = \frac{\mathcal{D}_a^{(0)} \langle \boldsymbol{\tau}_{ss} \rangle_h}{\mathcal{D}t} - (1+a) \operatorname{sym} \left( \left\langle \frac{\partial \boldsymbol{v}_s^{(1)}}{\partial x_3^{(1)}} \otimes \boldsymbol{\tau}_{3s} \right\rangle_h \right) \\ - \left[ \operatorname{sym} \left( \langle \{W_s(\boldsymbol{v}_s^{(1)}) + aD_s(\boldsymbol{v}_s^{(1)})\} \cdot \boldsymbol{\tau}_{ss} \rangle_h + (1+a) \left\langle \frac{\partial \boldsymbol{v}_s^{(2)}}{\partial x_3^{(1)}} \otimes \boldsymbol{\tau}_{3s} \right\rangle_h \right. \right. \\ \left. \left. + \frac{(1-a)h^{(1)}}{2} \nabla_s \operatorname{div}_s \boldsymbol{v}_s^{(0)} \otimes \langle \boldsymbol{\tau}_{3s} \rangle_h \right) \right] + \mathbf{O}(\epsilon^2), \quad (31a) $$

$$ \left\langle \left\langle \frac{\mathcal{D}_a \boldsymbol{\tau}}{\mathcal{D}t} \right\rangle_h \right\rangle_{s3} = \frac{\mathcal{D}_a^{(0)} \langle \boldsymbol{\tau}_{s3} \rangle_h}{\mathcal{D}t} + \frac{1-a}{2} \left\langle \frac{\partial \boldsymbol{v}_s^{(1)}}{\partial x_3^{(1)}} \cdot \boldsymbol{\tau}_{ss} \right\rangle_h - \frac{1+a}{2} \left\langle \frac{\partial \boldsymbol{v}_s^{(1)}}{\partial x_3^{(1)}} \tau_{33} \right\rangle_h + a \operatorname{div}_s (\boldsymbol{v}_s^{(0)}) \langle \boldsymbol{\tau}_{s3} \rangle_h \\ + \left\{ - \left\langle \left[ W_s(\boldsymbol{v}_s^{(1)}) + a D_s(\boldsymbol{v}_s^{(1)}) \right] \cdot \boldsymbol{\tau}_{s3} \right\rangle_h + \frac{h^{(1)}}{4} \left[ (1+a) \langle \boldsymbol{\tau}_{ss} \rangle_h - (1-a) \langle \boldsymbol{\tau}_{33} \rangle_h I_{ss} \right] \nabla_s \operatorname{div}_s \boldsymbol{v}_s^{(0)} \right. \\ \left. + \frac{1-a}{2} \left\langle \boldsymbol{\tau}_{ss} \cdot \frac{\partial \boldsymbol{v}_s^{(2)}}{\partial x_3^{(1)}} \right\rangle_h - \frac{1+a}{2} \left\langle \frac{\partial \boldsymbol{v}_s^{(2)}}{\partial x_3^{(1)}} \tau_{33} \right\rangle_h + a \langle \operatorname{div}_s (\boldsymbol{v}_s^{(1)}) \boldsymbol{\tau}_{3s} \rangle_h \right\} \epsilon + \mathbf{O}(\epsilon^2), \quad (31b) $$

$$ \left\langle \left\langle \frac{\mathcal{D}_a \boldsymbol{\tau}}{\mathcal{D}t} \right\rangle_h \right\rangle_{33} = \frac{D_s^{(0)} \langle \boldsymbol{\tau}_{33} \rangle_h}{Dt} + (1-a) \left\langle \frac{\partial \boldsymbol{v}_s^{(1)}}{\partial x_3^{(1)}} \cdot \boldsymbol{\tau}_{s3} \right\rangle_h + 2a \operatorname{div}_s (\boldsymbol{v}_s^{(0)}) \langle \boldsymbol{\tau}_{33} \rangle_h \\ + \left[ \left\langle \left( \frac{(1+a)h^{(1)}}{2} \nabla_s \operatorname{div}_s \boldsymbol{v}_s^{(0)} + (1-a) \frac{\partial \boldsymbol{v}_s^{(2)}}{\partial x_3^{(1)}} \right) \cdot \boldsymbol{\tau}_{s3} \right\rangle_h + 2a \langle \operatorname{div}_s (\boldsymbol{v}_s^{(1)}) \boldsymbol{\tau}_{33} \rangle_h \right] \epsilon + O(\epsilon^2), \quad (31c) $$

**Proof.** By noticing that we just need to develop the expression of the tensor product $(W(\boldsymbol{v}) + aD(\boldsymbol{v})) \cdot \boldsymbol{\tau}$ and to use the decomposition

$$ 2 \operatorname{sym}(\boldsymbol{\xi} \cdot \boldsymbol{\zeta}) = \left( \begin{array}{l} 2 \operatorname{sym}(\boldsymbol{\xi}_{ss} \cdot \boldsymbol{\zeta}_{ss} + \boldsymbol{\xi}_{s3} \otimes \boldsymbol{\zeta}_{3s}) \quad \boldsymbol{\xi}_{ss} \cdot \boldsymbol{\zeta}_{s3} + \boldsymbol{\xi}_{s3} \cdot \boldsymbol{\zeta}_{ss} + \boldsymbol{\xi}_{s3} \boldsymbol{\zeta}_{33} + \boldsymbol{\xi}_{33} \boldsymbol{\zeta}_{3s} \\ 2(\boldsymbol{\xi}_{3s} \cdot \boldsymbol{\zeta}_{s3} + \boldsymbol{\xi}_{33} \boldsymbol{\zeta}_{33}) \end{array} \right), $$

we can return to the proof of the previous theorem and conclude. □

<figure><img src="image_4.png" /></figure>
## 4. Examples

### 4.1. From Navier-Stokes equations to viscous shallow water equations

#### 4.1.1. Problem statement

The constitutive equation of an incompressible Newtonian fluid write:

$$ \sigma = 2\eta D(v) - pI, \qquad (32) $$

where $\sigma$ is the symmetric Cauchy stress tensor and $p$ a contribution to the pressure. The constant material parameter is the viscosity $\eta > 0$. These constitutive equations are coupled with standard conservation of mass (33a) and momentum (33b) equations. Assuming that the density $\rho$ is constant, i.e. the fluid is incompressible, the conservation equations write:

$$ \operatorname{div} v = 0, \qquad (33a) $$

$$ \rho \frac{Dv}{Dt} - \operatorname{div} \sigma = \rho g, \qquad (33b) $$

where $g = -g\hat{u}$ is the constant gravity vector, with $\hat{u}$ some given unit vector. The previous set of equations is closed by initial and boundary conditions. For the latter, we use, on the free surface, the conditions (4) with no surface tension ($\gamma = 0$), and, on the substrate, the conditions (6) with no friction ($k = 0$).

#### 4.1.2. Asymptotic analysis

With $\Sigma := \eta U/L$ as characteristic stress, the dimensionless problem writes:

**Problem 18.** Find $v$ defined in $\Omega$ such that

$$ \operatorname{Re} \frac{Dv}{Dt} - \operatorname{div} \sigma = - \frac{\operatorname{Re}}{Fr^2} \hat{u}, \qquad (34a) $$

$$ \operatorname{div} v = 0, \qquad (34b) $$

with $\sigma = 2D(v) - pI$, where $\operatorname{Re} := \rho UL/\eta$ is the Reynolds number and $\operatorname{Fr} := U/\sqrt{gL}$ is the Froude number.

The same notation is used for dimensioned and dimensionless field for convenience.

We would like to use the framework developed in Section 3.2 to get a shallow formulation of the above Navier-Stokes problem. To do so, $\mathcal{A}(u)$ is identified with $\operatorname{Re} Dv_i/Dt + (\operatorname{Re}/\operatorname{Fr}^2)\hat{u}_i$ and $j$ with $(\sigma_{ij})_{1\le j \le 3}$, for $i$ successively taken in $\{1,2,3\}$. Then, the condition of motion by slices is satisfied, up to order 1 at least, for each component of the velocity, thanks to Corollary 8. As a consequence, the rate of deformation expands with respect to $\epsilon$ as

$$ 2D(v) = \begin{pmatrix} 2D_s(v_s^{(0)}) & \mathbf{0} \\ \mathbf{0} & -2\operatorname{div}_s v_s^{(0)} \end{pmatrix} + \begin{pmatrix} 2D_s(v_s^{(1)}) & \frac{\partial v_s^{(2)}}{\partial x_3^{(1)}} - \nabla_s \operatorname{div}_s v_s^{(0)} x_3^{(1)} \\ \frac{\partial v_s^{(2)}}{\partial x_3^{(1)}} - \nabla_s \operatorname{div}_s v_s^{(0)\top} x_3^{(1)} & -2\operatorname{div}_s v_s^{(1)} \end{pmatrix} + \epsilon + \mathbf{O}(\epsilon^2), \quad (35) $$

where we used Proposition 12.

Following the literature, we make the fundamental assumption that $\operatorname{Re}/\operatorname{Fr}^2 = O(1/\epsilon)$. Therefore, due to the form of the deformation rate tensor (35), applying the Corollary 10 leads to

$$ p^{(-1)} = 0, \langle p^{(0)} \rangle_h = -2\operatorname{div}_s(v_s^{(0)}) + \frac{\hat{u}_3}{2} \frac{\epsilon \operatorname{Re}}{\operatorname{Fr}^2} h^{(1)}. \qquad (36) $$

We can now apply Theorem 7, Proposition 14 and Corollary 13 to get a shallow formulation of the above Navier-Stokes problem.
**Problem 19.** Find $\mathbf{v}_s^{(0)}$ and $h^{(1)}$ defined in $\Omega$ such that

$$ \frac{\partial h^{(1)}}{\partial t} + \operatorname{div}_s (h^{(1)} \mathbf{v}^{(0)}) = 0, \qquad (37a) $$

$$ \operatorname{Re} h^{(1)} \frac{D_s^{(0)} \mathbf{v}_s^{(0)}}{Dt} - \operatorname{div}_s (h^{(1)} \langle \boldsymbol{\sigma}_{ss}^{(0)} \rangle_h) = - \frac{\epsilon \operatorname{Re}}{\operatorname{Fr}^2} h^{(1)} \hat{\mathbf{u}}_s, \qquad (37b) $$

$$ \langle \boldsymbol{\sigma}_{ss}^{(0)} \rangle_h = 2D_s(\mathbf{v}_s^{(0)}) + 2\operatorname{div}_s(\mathbf{v}_s^{(0)}) \mathbf{I}_{ss} - \frac{\hat{\mathbf{u}}_3}{2} \frac{\epsilon \operatorname{Re}}{\operatorname{Fr}^2} h^{(1)} \mathbf{I}_{ss}. \qquad (37c) $$

This reduced model was first obtained in 2006 by Marche [10] after substantial formal calculations. In comparison, our approach requires less than one page. Solving Problem 19 instead of Problem 18 leads to an error of order $O(\epsilon)$ on velocity and height (cf. Theorem 7 or end of Section 3.3.1).

## 4.2. Maxwell viscoelastic model

We use the constitutive equations of a Maxwell viscoelastic fluid (1) along with initial and boundary conditions. For the latter, we use, on the free surface, the conditions (4) with no surface tension ($\gamma = 0$), and, on the substrate, the conditions (6) with no friction ($k = 0$).

With $\Sigma := \eta U/L$ again as characteristic stress, the dimensionless problem writes:

**Problem 20.** Find $\mathbf{v}$ and $\boldsymbol{\tau}$ defined in $\Omega$ such that

$$ \operatorname{Re} \frac{D\mathbf{v}}{Dt} - \operatorname{div}\boldsymbol{\sigma} = - \frac{\operatorname{Re}}{\operatorname{Fr}^2} \hat{\mathbf{u}}, \qquad (38a) $$

$$ \operatorname{div} \mathbf{v} = 0, \qquad (38b) $$

$$ \operatorname{We} \frac{\mathcal{D}_1 \boldsymbol{\tau}}{\mathcal{D} t} + \boldsymbol{\tau} = 2D(\mathbf{v}), \qquad (38c) $$

with $\boldsymbol{\sigma} = \boldsymbol{\tau} - p\mathbf{I}$, where We := $\lambda V/L$ is the Weissenberg number.

The same notation is used for dimensioned and dimensionless field for convenience.

### 4.2.1. Sufficient conditions for the motion by slices

With the help of our formalism, we can derive a new sufficient condition for the motion by slices in the context of the Maxwell model:

**Corollary 21 (motion by slices for the Maxwell model).** Let $\boldsymbol{\tau}$ be the elastic stress defined in the Maxwell model (1), and assume that it expands as $\boldsymbol{\tau} = \boldsymbol{\tau}^{(0)} + \boldsymbol{\tau}^{(1)}\epsilon + \mathbf{O}(\epsilon^2)$ and that the Weissenberg number We is of order $O(1)$, then

$$ \frac{\partial \mathbf{u}_s^{(0)}}{\partial x_3^{(1)}} = \frac{\partial \mathbf{u}_s^{(1)}}{\partial x_3^{(1)}} = \mathbf{0}. $$

**Proof.** Let us write the components s3 and 33 of the equation (38c):

$$ \operatorname{We} \left( \frac{D \boldsymbol{\tau}_{s3}}{Dt} - \tau_{33} \frac{\partial \mathbf{v}_s}{\partial x_3} - \left( \nabla_s \mathbf{v}_s + \frac{\partial \nu_3}{\partial x_3} \mathbf{I} \right) \cdot \boldsymbol{\tau}_{s3} - \boldsymbol{\tau}_{ss} \cdot \nabla_s \nu_3 \right) + \boldsymbol{\tau}_{s3} = \frac{\partial \mathbf{v}_s}{\partial x_3} + \nabla_s \nu_3, \qquad (39a) $$

$$ \operatorname{We} \left( \frac{D \tau_{33}}{Dt} - 2 \nabla_s \nu_3 \cdot \boldsymbol{\tau}_{s3} - 2 \frac{\partial \nu_3}{\partial x_3} \tau_{33} \right) + \tau_{33} = 2 \frac{\partial \nu_3}{\partial x_3}. \qquad (39b) $$

By assumption, we know the only term of order $O(1/\epsilon)$ in (39a) satisfies

$$ - \operatorname{We} \tau_{33}^{(0)} \frac{\partial \mathbf{v}_s^{(0)}}{\partial x_3^{(1)}} = \frac{\partial \mathbf{v}_s^{(0)}}{\partial x_3^{(1)}}. \qquad (39c) $$
Now, We = 0 leads to the result. Assume the opposite and that $\tau_{33}^{(0)} = -1/We$. Then, from equation (39b), we have

$$ \tau_{33}^{(0)} = -\frac{1}{We} = 0 \qquad (40) $$

since $-2We\tau_{33}^{(0)} = 2$ and $2\nabla_s v_3 \cdot \tau_{s3}$ is at most of order $O(\epsilon)$. This leads to a contradiction: $\tau_{33}^{(0)}$ cannot be uniformly constant equal to $-1/We$. We finally conclude from (39c) that $\partial v_s^{(0)} / \partial x_3^{(1)} = 0$. Note that if $\tau_{33}^{(0)} = 0$ is assumed, then the proof is still valid.

Hypothesis $\mathcal{H}^{(-1)}$ is satisfied then, from Lemma 5, there holds $\boldsymbol{\sigma}_{s3}^{(0)} = \mathbf{0}$, which leads to $\boldsymbol{\tau}_{s3}^{(0)} = \mathbf{0}$, and then, from (39a)

$$ (\mathrm{We}\tau_{33}^{(0)} + 1) \frac{\partial \boldsymbol{v}_s^{(1)}}{\partial x_3^{(1)}} = \mathbf{0}, \qquad (41) $$

which leads to the second result, since $\tau_{33}^{(0)}$ cannot be uniformly constant equal to $-1/We$. $\square$

### 4.2.2. Asymptotic analysis

In what follows, we assume that $\boldsymbol{\tau}$ asymptotic expansion is of the form $\boldsymbol{\tau} = \boldsymbol{\tau}^{(0)} + \boldsymbol{\tau}^{(1)}\epsilon + \boldsymbol{\tau}^{(2)}\epsilon^2 + \mathbf{O}(\epsilon^3)$, whereas $\mathrm{Re}/\mathrm{Fr}^2$ is of order $O(1/\epsilon)$, as in the previous example. Corollary 21 ensures that the property of motion by slices is satisfied.

**Momentum equation – Asymptotic.** The asymptotic analysis of the momentum equation has already been done in the previous example. Similar to what was done to deduce the pressure in the previous example, the zeroth order of the averaged pressure is given by $\langle p^{(0)} \rangle_h = \langle \tau_{33}^{(0)} \rangle_h + \hat{u}_3 \epsilon \mathrm{Re}/(2\mathrm{Fr}^2) h^{(1)}$, whence

$$ \langle \boldsymbol{\sigma}_{ss}^{(0)} \rangle_h = \langle \boldsymbol{\tau}_{ss}^{(0)} \rangle_h - \langle \boldsymbol{\tau}_{33}^{(0)} \rangle_h \mathbf{I}_{ss} - \frac{\hat{u}_3}{2} \frac{\epsilon \mathrm{Re}}{\mathrm{Fr}^2} h^{(1)} \mathbf{I}_{ss}. \qquad (42) $$

In addition, assumption $\mathcal{H}^{(-1)}$ (see Lemma 5) is satisfied for both planar velocity components, then

$$ \boldsymbol{\sigma}_{s3}^{(0)} = \boldsymbol{\tau}_{s3}^{(0)} = \mathbf{0}. \qquad (43) $$

**Maxwell equation – Asymptotic.** Taking into account the very last result (43) and applying Theorem 17 and Proposition 12 lead to the shallow formulation of the Maxwell equation (1):

$$ \mathrm{We} \frac{\mathcal{D}_1^{(0)} \langle \boldsymbol{\tau}_{ss}^{(0)} \rangle_h}{\mathcal{D} t} + \langle \boldsymbol{\tau}_{ss}^{(0)} \rangle_h = 2D_s (\boldsymbol{v}_s^{(0)}), \qquad (44a) $$

$$ \mathrm{We} \left( \frac{D_s^{(0)} \langle \tau_{33}^{(0)} \rangle_h}{Dt} + 2 \mathrm{div}_s (\boldsymbol{v}_s^{(0)}) \langle \tau_{33}^{(0)} \rangle_h \right) + \langle \tau_{33}^{(0)} \rangle_h = -2 \mathrm{div}_s (\boldsymbol{v}_s^{(0)}). \qquad (44b) $$

**Remark 22.** From a microscopic perspective, it is known [3, p. 90] that the elastic stress can be expressed in terms of the conformation tensor $\boldsymbol{c} \propto \langle \boldsymbol{\ell} \otimes \boldsymbol{\ell} \rangle$ defined as an average over all possible directions $\boldsymbol{\ell}$ of molecules constituting the considered fluid:

$$ \boldsymbol{\tau} = \frac{\eta}{\lambda} (\boldsymbol{c} - \mathbf{I}). \qquad (45) $$

On the reasonable assumption that $\boldsymbol{\ell}$ expands with respect to $\epsilon$ as $\boldsymbol{\ell} = \begin{pmatrix} \mathbf{o}^{(1)} \\ O(\epsilon) \end{pmatrix}$, we deduce a possible expansion for the conformation tensor as $\boldsymbol{c} = \begin{pmatrix} \mathbf{o}^{(1)} & \mathbf{o}(\epsilon) \\ \mathbf{o}(\epsilon) & O(\epsilon^2) \end{pmatrix}$. Thus, at zeroth order, the component 33 of the elastic stress reads $\tau_{33}^{(0)} = -1/We$, but it would lead to the contradiction $-1/We = 0$, according to equation (44b). If this expansion does not work for $\boldsymbol{c}$, does it at least work for $\boldsymbol{\tau}$? Assume $\tau_{33}^{(0)} = 0$, then $\mathrm{div}_s(\boldsymbol{v}_s^{(0)}) = 0$ from equation (44b). If this does not lead to a contradiction, it implies a strong physical constraint, which does not hold in most situations. We therefore understand that the order with respect to $\epsilon$ of a given quantity (here $\boldsymbol{\tau}$) must be componentwise at least of the order of the second member (here $D(\boldsymbol{v})$).
**Summary.** The shallow formulation of the Maxwell problem writes

**Problem 23.** Find $\boldsymbol{v}_s^{(0)}$, $h^{(1)}$, $\langle \boldsymbol{\tau}_{ss}^{(0)} \rangle_h$ and $\langle \tau_{33}^{(0)} \rangle_h$ defined in $\Omega$ such that

$$ \frac{\partial h^{(1)}}{\partial t} + \mathrm{div}_s(h^{(1)} \boldsymbol{v}^{(0)}) = 0, \qquad (46a) $$

$$ \mathrm{Re} h^{(1)} \frac{D_s^{(0)} \boldsymbol{v}_s^{(0)}}{Dt} - \mathrm{div}_s(h^{(1)} \langle \boldsymbol{\sigma}_{ss}^{(0)} \rangle_h) = - \frac{\epsilon \mathrm{Re}}{\mathrm{Fr}^2} h^{(1)} \hat{\boldsymbol{u}}_s, \qquad (46b) $$

$$ \langle \boldsymbol{\sigma}_{ss}^{(0)} \rangle_h = \langle \boldsymbol{\tau}_{ss}^{(0)} \rangle_h - \langle \tau_{33}^{(0)} \rangle_h \boldsymbol{I}_{ss} - \frac{\hat{u}_3}{2} \frac{\epsilon \mathrm{Re}}{\mathrm{Fr}^2} h^{(1)} \boldsymbol{I}_{ss}, \qquad (46c) $$

$$ \mathrm{We} \frac{\mathcal{D}_1^{(0)} \langle \boldsymbol{\tau}_{ss}^{(0)} \rangle_h}{\mathcal{D} t} + \langle \boldsymbol{\tau}_{ss}^{(0)} \rangle_h = 2 D_s (\boldsymbol{v}_s^{(0)}), \qquad (46d) $$

$$ \mathrm{We} \left( \frac{D_s^{(0)} \langle \tau_{33}^{(0)} \rangle_h}{Dt} + 2 \mathrm{div}_s(\boldsymbol{v}_s^{(0)}) \langle \tau_{33}^{(0)} \rangle_h \right) + \langle \tau_{33}^{(0)} \rangle_h = -2 \mathrm{div}_s(\boldsymbol{v}_s^{(0)}). \qquad (46e) $$

It is consistent with what has been derived in [4]. When $\mathrm{We} = 0$, it reduces to the same asymptotic (37) as the Navier–Stokes equation. Solving Problem 23 instead of Problem 20 leads to an error of order $O(\epsilon)$ on velocity, height and all other averaged quantities (cf. Theorem 7 or end of Section 3.3.1).

## 5. Conclusion

A new framework for asymptotic analysis has been presented in this paper. We have proposed a set of tools for performing thin-layer approximations of free-surface fluid models with flat topography and no-penetration boundary condition. The motion by slices is no more an arbitrary assumption, and minimal assumptions for obtaining it have been pointed out. Moreover, for first time, to the best of our knowledge, the first-order expansion of the Gordon–Schowalter derivatives has been developed. The proposed framework makes it possible to clarify the necessary assumptions and to simplify the asymptotic analysis of a wide variety of partial differential equations from continuum physics, whether they are scalar, vector or tensor, with objective derivative and/or diffusive, reactive or advective terms. This framework has been illustrated by two applications: the Navier–Stokes and the Maxwell models. In perspective, a generalisation of this work would be to consider here an arbitrary topography, following the work made in [6] for the Saint-Venant system.

## Appendix A. Leibniz formulas

**Proposition 24 (Leibniz's Integral Rule).** Let $\mathcal{X}$ and $I$ be two intervals of $\mathbb{R}$ and $f: (x, t) \in \mathcal{X} \times I \to f(x, t) \in \mathbb{R}^N$, $N \in \mathbb{N}\{0\}$, a sufficiently smooth function. Let also $a$ and $b$ two differentiable functions defined from $\mathcal{X}$ to $I$. Then the parameterized integral function $F$, defined in $\mathcal{X}$ by

$$ F(x) = \int_{a(x)}^{b(x)} f(x, t) dt, $$

is differentiable and

$$ F'(x) = \int_{a(x)}^{b(x)} \frac{\partial f}{\partial x} dt + f(x, b(x)) b'(x) - f(x, a(x)) a'(x). $$

**Corollary 25 (Leibniz's Integral Rule for depth average).** Let $\mathcal{L}$ be a first order differential linear operator, depending only on partial derivatives $\partial_t$, $\partial_{x_1}$ and $\partial_{x_2}$. For any sufficiently smooth scalar field $\varphi$, the following formula is satisfied:

$$ \mathcal{L}(h\langle\varphi\rangle_h) = h\langle\mathcal{L}\varphi\rangle_h + \varphi(x_3 = h)\mathcal{L}h. \qquad (47a) $$
**Corollary 26.** For any sufficiently smooth vector field **u**, the following formulas are satisfied:

$$ \mathrm{div}_s(h\langle\mathbf{u}\rangle_h) = h\langle\mathrm{div}_s\mathbf{u}\rangle_h + \mathbf{u}_s(x_3 = h) \cdot \nabla_s h, \quad (47b) $$

$$ \Delta_s(h\langle\mathbf{u}\rangle_h) = h\langle\Delta_s\mathbf{u}\rangle_h + 2\nabla_s\mathbf{u}_s(x_3 = h) \cdot \nabla_s h + \mathbf{u}_s(x_3 = h)\Delta_s h, \quad (47c) $$

$$ h\langle\Delta_s\mathbf{u}\rangle_h = \mathrm{div}_s(h\langle\nabla_s\mathbf{u}\rangle_h) - \nabla_s\mathbf{u}_s(x_3 = h) \cdot \nabla_s h. \quad (47d) $$

We close this section by giving a last relation that will be useful:

$$ \langle x_3 \partial_{x_3} \varphi \rangle_h + \langle \varphi \rangle_h = \varphi(x_3 = h). \quad (47e) $$

## Acknowledgements

This work has taken place in the context of an interdisciplinary collaboration with two physicists, Hélène Delanoë-Ayari (Institut Lumière Matière, Univ. Lyon 1) and François Graner (Laboratoire Matière et Systèmes Complexes, Univ. Paris 7), on the modelling of biological tissues. We warmly thank them for fruitful discussions. This work is supported by the CNRS interdisciplinary mission (MITI).

## References

[1] N. J. Balmforth, R. V. Craster, "A consistent thin-layer theory for Bingham plastics", *J. Non-Newton. Fluid Mech.* **84** (1999), no. 1, p. 65-81.

[2] N. Bernabeu, P. Saramito, C. Smutek, "Numerical modeling of non-Newtonian viscoplastic flows. II: Viscoplastic fluids and general tridimensional topographies", *Int. J. Numer. Anal. Model.* **11** (2014), no. 1, p. 213-228.

[3] R. B. Bird, C. F. Curtiss, R. C. Armstrong, O. Hassager, *Dynamics of polymeric liquids. Volume 2: Kinetic Theory*, 2nd ed., Wiley-Interscience, 1987.

[4] F. Bouchut, S. Boyaval, "A new model for shallow viscoelastic fluids", *Math. Models Methods Appl. Sci.* **23** (2013), no. 8, p. 1479-1526.

[5] ——, "Unified derivation of thin-layer reduced models for shallow free-surface gravity flows of viscous fluids", *Eur. J. Mech. B Fluids* **55** (2016), p. 116-131.

[6] F. Bouchut, M. Westdickenberg, "Gravity driven shallow water models for arbitrary topography", *Commun. Math. Sci.* **2** (2004), no. 3, p. 359-389.

[7] E. D. Fernández-Nieto, P. Noble, J.-P. Vila, "Shallow water equations for non-Newtonian fluids", *J. Non-Newton. Fluid Mech.* **165** (2010), no. 13, p. 712-732.

[8] J.-F. Gerbeau, B. Perthame, "Derivation of Viscous Saint-Venant System for Laminar Shallow Water; Numerical Validation", *Research Report RR-4084*, INRIA, 2000, Projet M3N.

[9] K. F. Liu, C. C. Mei, "Approximate equations for the slow spreading of a thin sheet of Bingham plastic fluid", *Phys. Fluids, A2* (1990), no. 1, p. 30-36.

[10] F. Marche, "Derivation of a new two-dimensional viscous shallow water model with varying topography, bottom friction and capillary effects", *Eur. J. Mech. B Fluids* **26** (2007), no. 1, p. 49-63.

[11] J. G. Oldroyd, "On the formulation of rheological equations of state", *Proc. R. Soc. Lond., Ser. A 200* (1950), p. 523-541.

[12] A. B. de Saint-Venant, "Théorie et équations générales du mouvement non permanent des eaux courantes", *C. R. Acad. Sci., Paris 17* (1871), no. 73, p. 147-154.