---
title: Backward Euler Stokes–Biot Splitting
url: https://www.emergentmind.com/topics/backward-euler-stokes-biot-splitting-method
type: topic
---

# Backward Euler Stokes–Biot Splitting

Searching arXiv for the cited Stokes–Biot backward-Euler splitting literature to ground the article in the relevant papers.
arXiv search query: "Stokes Biot splitting backward Euler Robin Robin fluid poroelastic interaction"
The backward Euler Stokes–Biot splitting method denotes a class of fully discrete time-stepping procedures in which a poroelastic Biot model, or a coupled Stokes–Biot fluid–poroelastic system, is advanced by backward Euler in time while the multiphysics coupling is decomposed into smaller subproblems. In one line of work, the Biot equations are reformulated into a generalized Stokes problem coupled with a diffusion problem, yielding a two-step finite element method with a discrete energy law and optimal error estimates [2204.03201]. In another line, transient Stokes flow and poroelasticity are decoupled across an interface by explicit time-lagging, Robin transmission operators, or interface Lagrange multipliers, so that the fluid and poroelastic subproblems can be solved separately at each time step [2603.21088], [2409.18910], [2601.14201]. Extensions also combine the backward Euler Stokes–Biot split with fixed-strain Biot decoupling for multilayered poroelastic structures interacting with Stokes flow [2507.10538].

## 1. Continuous models and the meaning of the split

In the reformulated Biot setting with secondary consolidation, the unknowns are the solid displacement $u:\Omega\times[0,T]\to\mathbb{R}^d$ and pore pressure $p:\Omega\times[0,T]\to\mathbb{R}$ on a polygonal domain $\Omega\subset\mathbb{R}^d$, $d=2,3$. The constitutive law uses
$$
\sigma(u)=2\mu \epsilon(u)+\lambda\,\mathrm{tr}(\epsilon(u))I,\qquad 
\epsilon(u)=\tfrac12(\nabla u+\nabla u^T),
$$
with parameters $\lambda^\ast,a_0,b_0,\mu,\lambda>0$, positive-definite $K$, and $\theta_f>0$. The key algebraic reformulation introduces
$$
q:=\nabla\cdot u,\qquad \Phi:=a_0 p+b_0 q,\qquad \Delta:=b_0 p-q-\lambda^\ast q_t,
$$
which converts the original fluid–solid coupled problem into a Stokes-like system for $(u,\Delta)$ coupled to a scalar diffusion problem in $\Phi$ [2204.03201]. In that formulation,
$$
-\mu \Delta u+\Delta =F,
$$
is coupled to an equation relating $\Delta$, $q$, $q_t$, and $\Phi$, and to a diffusion equation for $\Phi_t$ involving $K$, $\theta_f$, and the reconstructed pore pressure.

In fluid–poroelastic interaction, the split has a different meaning. The computational domain is partitioned into a fluid region $\Omega_f$ governed by the transient Stokes equations and a poroelastic region $\Omega_p$ governed by the Biot equations, with interface conditions enforcing mass conservation, stress balance, and the Beavers–Joseph–Saffman law [2409.18910], [2601.14201]. In the multilayered variant, a thick Biot layer $\Omega_e$ is coupled to a thin poroelastic plate and to a Stokes fluid domain $\Omega_f$; the interface conditions include continuity of normal velocity, continuity of displacement, balance of normal forces, the Beavers–Joseph–Saffman slip law, and continuity of poro-pressures [2507.10538].

Across these formulations, the phrase “Stokes–Biot splitting” therefore refers either to an internal decomposition of the Biot model into generalized Stokes and diffusion components, or to a partitioned treatment of a coupled Stokes/Biot system. The common denominator is backward Euler time integration combined with decoupling of the multiphysics operators.

## 2. Reformulated Biot method: generalized Stokes step and diffusion step

For the Biot model with secondary consolidation, the finite element discretization is built on
$$
V:=H^1(\Omega)^d,\qquad Q:=L^2(\Omega),\qquad Q_0:=\{\phi\in Q;\int_\Omega \phi=0\}.
$$
On a quasi-uniform triangulation $\mathcal{T}_h$, the method uses the Taylor–Hood pair for the generalized Stokes variables and piecewise linear elements for the diffusion variable:
$$
X_h:=\{v_h\in C^0; \ v_h|_K\in P_2(K)^d\},\qquad
M_h:=\{\phi_h\in C^0; \ \phi_h|_K\in P_1(K)\},
$$
with $W_h:=M_h$, $M_{h,0}=M_h\cap Q_0$, and
$$
V_h:=\{v_h\in X_h;(v_h,r)=0 \ \forall \text{ rigid motions } r\}.
$$
The pair $(X_h,M_h)$ satisfies the discrete inf-sup condition for $\nabla\cdot v_h$ [2204.03201].

With time grid $0=t_0<t_1<\cdots<t_N=T$, $\Delta t=T/N$, one step of the fully discrete backward-Euler splitting proceeds in two stages. The first is a generalized Stokes subproblem: given $(u_h^n,\Delta_h^n,\Phi_h^n)$, solve for $(u_h^{n+1},\Delta_h^{n+1})\in V_h\times M_h$ so that
$$
\mu(\epsilon(u_h^{n+1}),\epsilon(v_h))-(\Delta_h^{n+1},\nabla\cdot v_h)
=(F,v_h)+\langle F_1,v_h\rangle,
$$
and
$$
\chi_3(\Delta_h^{n+1},\psi_h)+(q_h^{n+1},\psi_h)+\lambda^\ast\chi_3(d_tu_h^{n+1},\psi_h)
=\chi_1(\Phi_h^n,\psi_h),
$$
where $d_tu_h^{n+1}=(u_h^{n+1}-u_h^n)/\Delta t$ and $q_h^{n+1}:=\nabla\cdot u_h^{n+1}$. The second is a diffusion subproblem: solve for $\Phi_h^{n+1}\in W_h$ from
$$
(d_t\Phi_h^{n+1},\eta_h)+\frac{1}{\theta_f}(K\nabla p_h^{n+1},\nabla \eta_h)
=(\phi,\eta_h)+\langle \phi_1,\eta_h\rangle,
$$
with reconstructed pressure
$$
p_h^{n+1}=\chi_1\Delta_h^{n+1}+\chi_2\Phi_h^{n+1}+\lambda^\ast\chi_1 d_t q_h^{n+1}.
$$
The right-hand side in the generalized Stokes equation uses $\Phi_h^n$; a symmetric variant may instead use $\Phi_h^{n+\frac12}$ [2204.03201].

Initialization is explicit: $u_h^0$ is the projection of $u_0$, $\Delta_h^0$ the projection of $b_0p_0-q_0$, and $\Phi_h^0=a_0p_0+b_0q_0$. Each time step then consists of solving the generalized Stokes problem, updating $q_h^{n+1}=\nabla\cdot u_h^{n+1}$, solving the scalar diffusion problem, and reconstructing $p_h^{n+1}$ algebraically [2204.03201].

## 3. Interface-partitioned Stokes–Biot schemes

For coupled fluid–poroelastic interaction, the backward-Euler split is usually imposed at the interface rather than through an algebraic reformulation of the Biot operator. The resulting methods differ mainly in how interface data are transferred between the Stokes and Biot solves.

| Formulation | Decoupling mechanism | Reported property |
|---|---|---|
| Explicit splitting [2603.21088] | old interface traces in both subproblems | uncoupled variational problems solved in parallel |
| Robin–Robin splitting [2409.18910] | auxiliary Robin data $\lambda^n=(\mu_n^n,\mu_\tau^n)$ | single decoupled Stokes and Biot solves per time step |
| Lagrange multiplier partitioned method [2601.14201] | Schur complement in interface pressure and multipliers | fluid and poroelastic subproblems solved independently at each time step |
| Multilayered split [2507.10538] | backward Euler Stokes–Biot split plus fixed-strain Biot splitting | three independent linear solves per time step |

In the explicit scheme, the fluid subproblem and the poroelastic subproblem are solved independently at each time step, using interface data from the previous time level. The Stokes step employs Nitsche-type penalty parameters $\gamma>0$ and $L>0$ to enforce tangential-slip and normal-flux/pressure continuity, while all right-hand interface data appear at time level $n$, making the coupling explicit [2603.21088]. The stated consequence is that the fluid and Biot subproblems can be solved in parallel.

In the Robin–Robin method, the interface information is represented by an auxiliary variable $\lambda^n=(\mu_n^n,\mu_\tau^n)\in\Lambda_h$, which approximates the Robin boundary data and is updated after one Stokes solve and one Biot solve. The fluid subproblem weakly imposes Robin conditions on normal and tangential components of the Stokes stress, and the poroelastic subproblem weakly imposes Robin conditions involving displacement, velocity, and pore pressure; the update of $\lambda^{n+1}$ transfers the latest stress–velocity mismatch to the next time step [2409.18910].

In the Lagrange multiplier partitioned method, the coupled weak formulation introduces three interface multipliers,
$$
\lambda_1:=(\sigma_f n_f)\cdot n_f,\qquad
\lambda_2:=(\sigma_f n_f)\cdot \tau,\qquad
\lambda_p:=K\nabla p_p\cdot n_p,
$$
to enforce mass balance, stress continuity, and the Beavers–Joseph–Saffman law. After backward-Euler discretization and finite element assembly, the monolithic system is reduced by block Gaussian elimination to a Schur complement problem for $(p_f^{n+1},\lambda^{n+1})$, after which the fluid velocity, pore pressure, and displacement are recovered by forward solves [2601.14201].

The multilayered formulation combines two distinct splittings. First, the Stokes problem is decoupled from the multilayered structure problem. Second, within the poroelastic structure, the fixed-strain Biot split decouples flow and mechanics. The resulting time step consists of a Darcy-type solve for the thick Biot domain and thin poroelastic plate, an elastodynamics solve for the thick solid and plate variables, and a final Stokes solve driven by the updated interface traction [2507.10538].

## 4. Discrete energy laws and stability mechanisms

A central feature of backward-Euler Stokes–Biot splitting is that stability is typically expressed through a discrete energy identity or energy inequality, but the form of that estimate depends strongly on the splitting strategy.

For the reformulated Biot method, the discrete energy law is exact for $\theta=0$ or $1$:
$$
J_{h,\theta}^{\ell+1}+S_{h,\theta}^{\ell+1}=J_{h,\theta}^0.
$$
Here
$$
J_{h,\theta}^{n+1}
=\tfrac12\Big[\mu\|\epsilon(u_h^{n+1})\|^2
+\chi_2\|\Phi_h^{n+\theta}\|^2
+\chi_3\|\lambda^\ast d_tu_h^{n+1}+\Delta_h^{n+1}\|^2
+\text{terms of }O(\Delta t)^2\Big],
$$
while $S_{h,\theta}^{\ell+1}$ is a cumulative dissipation containing terms such as
$$
\Delta t\sum_{n=0}^{\ell}\Big[\lambda^\ast\|d_tu_h^{n+1}\|^2
+\frac1{\theta_f}(K\nabla p_h^{n+1},\nabla p_h^{n+1})+\cdots\Big].
$$
For $\theta=0$, the method yields a discrete non-increase under the mesh–CFL condition $\Delta t=O(h^2)$ [2204.03201].

For the explicit Stokes–Biot interface split, the analysis introduces a total error energy $X_n^2$ and a dissipative quantity $Y_n^2$, then proves the exact energy–dissipation balance
$$
X_n^2-X_0^2+\sum_{i=1}^n Y_i^2+\sum_{i=1}^n Z_i
=\sum_{i=1}^n\big[F_u^i(e_u^{h,i})+F_\xi^i(e_\xi^{h,i})+F_\phi^i(e_\phi^{h,i})\big].
$$
After bounding the mixed interface terms $Z_i$ and the consistency terms $F_\bullet^i$, one obtains a discrete Gronwall inequality and hence unconditional stability in a combined energy-dissipation norm [2603.21088].

The Robin–Robin scheme has an unconditional stability estimate of the form
$$
E^N+\Delta t\sum_{n=1}^N D^n+\sum_{n=1}^N S^n\le E^0,
$$
for zero volume forces, independently of $\Delta t$ and $h$. The corresponding energy includes fluid kinetic energy, poroelastic kinetic and elastic energy, pore-pressure energy, and an interface contribution involving the auxiliary Robin data $\lambda^n$ [2409.18910].

By contrast, the multilayered split is proved conditionally stable. With zero external forcing and a CFL-type condition on $\Delta t$ involving mesh size $h$, thickness $H$, and physical parameters, the scheme satisfies
$$
\mathcal{E}^{n+1}+\epsilon(\|\eta^{n+1}-\eta^n\|^2+\cdots)+2\Delta t\,\mathcal{D}^{n+1}\le \mathcal{E}^n.
$$
Under the further parameter constraints $\alpha^2<c_0\lambda_e$ and $(\alpha^p)^2<12c_0^pD$, the paper states a milder CFL condition $\Delta t\lesssim \min\{h,\ldots\}$ and unconditional stability in the structural coupling parameters [2507.10538].

These results distinguish a common source of confusion. Backward Euler does not by itself determine whether a Stokes–Biot partitioned method is unconditional: unconditional estimates are proved for some interface treatments, whereas other decompositions require CFL-type restrictions.

## 5. Error estimates and convergence orders

For the reformulated Biot method with $P_2$–$P_1$–$P_1$ elements, the a priori estimate states that, under sufficient regularity, for $\theta=0$ with $\Delta t=O(h^2)$ or for $\theta=1$ unconditionally,
$$
\max_{0\le k\le n}\Big[
\|\epsilon(u(t_{k+1})-u_h^{k+1})\|
+\|\Phi(t_{k+1})-\Phi_h^{k+1}\|
+\|\Delta(t_{k+1})-\Delta_h^{k+1}\|
\Big]
+\Big(\Delta t\sum_{k=0}^n \|\nabla(p-p_h)\|^2\Big)^{1/2}
\le C(\Delta t+h^2).
$$
The same source states that the spatial discretization delivers second-order convergence in $L^2$ for pressures and third-order in $L^2$ for displacement, while backward Euler gives first-order convergence in time [2204.03201]. Numerical tests on $\Omega=(0,1)^2$ report
- $L^2(\Omega)$-error$(u_h)\simeq O(h^3)$ and $H^1$-error$(u_h)\simeq O(h^2)$,
- $L^2$-error$(p_h)\simeq O(h^2)$ and $H^1$-error$(p_h)\simeq O(h)$,
- first-order in $\Delta t$,
with representative rates $3.04$, $2.09$, $3.02$, $2.03$, $2.24$, and $1.03$ in the reported manufactured-solution tables [2204.03201].

For the explicit interface split, the main estimate is
$$
\max_{0\le n\le N} X_n+\Big(\sum_{n=1}^N Y_n^2\Big)^{1/2}
\le C\big(\Delta t+h^{k+r/2}\big),
$$
where $k$ is the polynomial degree used in the finite element spaces and $r\in(0,1]$ is the extra regularity exponent in the $L^2$-dual estimates, with $r=1$ on convex domains. The stated interpretation is first-order accuracy in time and optimal spatial convergence rates determined by the approximation spaces [2603.21088].

For the Lagrange multiplier partitioned method, the convergence theorem gives
$$
\max_{0\le m\le M}\Big\{
\|u(t_m)-u_h^m\|_{H^1(\Omega_f)}^2
+\|n(t_m)-n_h^m\|_{H^1(\Omega_p)}^2
+\|p_f(t_m)-p_{f,h}^m\|_{L^2(\Omega_f)}^2
+\|p_p(t_m)-p_{p,h}^m\|_{L^2(\Omega_p)}^2
+\text{interface-multiplier norms}
\Big\}^{1/2}
\le C(\Delta t+h^k),
$$
with $k$ dictated by the polynomial degrees [2601.14201].

For the Robin–Robin split, the time-discretization error satisfies
$$
\max_{1\le n\le N}\Big(
\|\eta(t_n)-\eta^n\|_e+\|p_p(t_n)-p_p^n\|_{L^2(\Omega_p)}
\Big)
+\Big(\Delta t\sum_{n=1}^N\|u_f(t_n)-u_f^n\|_f^2\Big)^{1/2}
+\cdots
\le C\sqrt{T}\,\Delta t,
$$
and the paper states the time discretization error as $\mathcal{O}(\sqrt{T}\Delta t)$ [2409.18910].

For the multilayered split, the numerical study reports spatial convergence $O(h^2)$ in $L^2$-norm for all variables when using $P1$ finite elements for poroelastic variables and $P2$–$P1$ Taylor–Hood for the fluid, together with temporal convergence $O(\Delta t)$ for all variables [2507.10538]. Since these estimates arise in different norms and for different interface formulations, they are complementary rather than directly comparable.

## 6. Variants, implementation, and interpretive context

Several implementation patterns recur throughout the literature. The reformulated Biot method leads to a two-step algorithm whose unknowns are $(u_h,\Delta_h,\Phi_h)$ and from which the pore pressure is reconstructed algebraically [2204.03201]. The explicit interface split leads to two uncoupled variational solves per time step, one in $\Omega_f$ and one in $\Omega_p$, with all cross-domain traces lagged to time level $n$ [2603.21088]. The Robin–Robin scheme similarly performs one Stokes solve, one Biot solve, and one interface update, while an iterative version repeats these substeps until convergence and is shown to converge to the fully implicit monolithic backward-Euler scheme [2409.18910].

The multilayered method is structurally more elaborate. Its Darcy-type substep computes the thick-domain pore pressure $p^{n+1}$, plate pressure $q^{n+1}$, and Darcy velocities $u_e^{n+1}$ and $u_p^{n+1}$ using fixed-strain data $\xi^n$, $v^n$, and interface fluid velocity $u^n$. Its structure substep advances $\eta^{n+1}$, $\xi^{n+1}$, $w^{n+1}$, $v^{n+1}$, and $\Lambda^{n+1}$ from the updated pressures and interface traction. Its fluid substep then solves the Stokes system for $(u^{n+1},\pi^{n+1})$ with tangential slip and normal traction prescribed by the new plate–poroelastic data [2507.10538]. The reported implementation uses FEniCS and standard direct or iterative solvers for each subproblem, together with an additional penalty term
$$
(\gamma_{\mathrm{pen}}\mu_f/d_h)\langle u\cdot n-v-u_p,\ u\cdot n-v-u_p\rangle_\Gamma
$$
to improve mass-conservation [2507.10538].

The Schur-complement formulation provides a different computational perspective. Instead of alternating explicit interface data, it keeps the backward-Euler monolithic discretization intact but eliminates bulk variables, leaving an interface system for $(p_f^{n+1},\lambda^{n+1})$. Matrix-vector products with the Schur complement are evaluated through sparse solves with time-invariant matrices $W_f$, $W_n$, and $W_p$, and the preconditioner
$$
M_{\rm pre}:=A_0+\Delta t\,A_1
$$
is reported to reduce the spectrum of the Schur complement to $O(1)$ and to yield iteration counts $O(1)$ independent of mesh size [2601.14201].

Taken together, these formulations show that “backward Euler Stokes–Biot splitting” is not a single canonical algorithm but a family of backward-Euler partitionings. The decisive design choices are the location of the split—inside the Biot operator, across the Stokes/Biot interface, or within a multilayer poroelastic substructure—and the interface transfer mechanism, which may be explicit, Robin-based, or multiplier-based. The cited analyses indicate that these choices determine whether the method is exactly energy balanced, unconditionally stable, conditionally stable under a CFL-type restriction, or optimized for Schur-complement solution and preconditioning rather than for explicit loose coupling.

Source: https://www.emergentmind.com/topics/backward-euler-stokes-biot-splitting-method