Backward Euler Stokes–Biot Splitting

• Backward Euler Stokes–Biot splitting is a family of partitioned algorithms that uses backward Euler time integration combined with a decomposition of coupled fluid and poroelastic problems.
• The method separates the original problem either by reformulating the Biot system into a generalized Stokes plus diffusion formulation or by decoupling fluid and structure through explicit time-lagging, Robin conditions, or Lagrange multipliers.
• Discrete energy laws and error estimates show that the schemes achieve first-order time convergence and optimal spatial convergence while ensuring stability under various CFL-type conditions or unconditionally.

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 (Ge et al., 2022). 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 (Wang et al., 22 Mar 2026, Dalal et al., 2024, Castro et al., 20 Jan 2026). Extensions also combine the backward Euler Stokes–Biot split with fixed-strain Biot decoupling for multilayered poroelastic structures interacting with Stokes flow (Scharf et al., 14 Jul 2025).

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 $p:\Omega\times[0,T]\to\mathbb{R}$0 (Ge et al., 2022). In that formulation,

$p:\Omega\times[0,T]\to\mathbb{R}$1

is coupled to an equation relating $p:\Omega\times[0,T]\to\mathbb{R}$2, $p:\Omega\times[0,T]\to\mathbb{R}$3, $p:\Omega\times[0,T]\to\mathbb{R}$4, and $p:\Omega\times[0,T]\to\mathbb{R}$5, and to a diffusion equation for $p:\Omega\times[0,T]\to\mathbb{R}$6 involving $p:\Omega\times[0,T]\to\mathbb{R}$7, $p:\Omega\times[0,T]\to\mathbb{R}$8, and the reconstructed pore pressure.

In fluid–poroelastic interaction, the split has a different meaning. The computational domain is partitioned into a fluid region $p:\Omega\times[0,T]\to\mathbb{R}$9 governed by the transient Stokes equations and a poroelastic region $\Omega\subset\mathbb{R}^d$0 governed by the Biot equations, with interface conditions enforcing mass conservation, stress balance, and the Beavers–Joseph–Saffman law (Dalal et al., 2024, Castro et al., 20 Jan 2026). In the multilayered variant, a thick Biot layer $\Omega\subset\mathbb{R}^d$1 is coupled to a thin poroelastic plate and to a Stokes fluid domain $\Omega\subset\mathbb{R}^d$2; 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 (Scharf et al., 14 Jul 2025).

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

$\Omega\subset\mathbb{R}^d$3

On a quasi-uniform triangulation $\Omega\subset\mathbb{R}^d$4, the method uses the Taylor–Hood pair for the generalized Stokes variables and piecewise linear elements for the diffusion variable:

$\Omega\subset\mathbb{R}^d$5

with $\Omega\subset\mathbb{R}^d$6, $\Omega\subset\mathbb{R}^d$7, and

$\Omega\subset\mathbb{R}^d$8

The pair $\Omega\subset\mathbb{R}^d$9 satisfies the discrete inf-sup condition for $d=2,3$0 (Ge et al., 2022).

With time grid $d=2,3$1, $d=2,3$2, one step of the fully discrete backward-Euler splitting proceeds in two stages. The first is a generalized Stokes subproblem: given $d=2,3$3, solve for $d=2,3$4 so that

$d=2,3$5

and

$d=2,3$6

where $d=2,3$7 and $d=2,3$8. The second is a diffusion subproblem: solve for $d=2,3$9 from

$$\sigma(u)=2\mu \epsilon(u)+\lambda\,\mathrm{tr}(\epsilon(u))I,\qquad \epsilon(u)=\tfrac12(\nabla u+\nabla u^T),$$0

with reconstructed pressure

$$\sigma(u)=2\mu \epsilon(u)+\lambda\,\mathrm{tr}(\epsilon(u))I,\qquad \epsilon(u)=\tfrac12(\nabla u+\nabla u^T),$$1

The right-hand side in the generalized Stokes equation uses $$\sigma(u)=2\mu \epsilon(u)+\lambda\,\mathrm{tr}(\epsilon(u))I,\qquad \epsilon(u)=\tfrac12(\nabla u+\nabla u^T),$$2; a symmetric variant may instead use $$\sigma(u)=2\mu \epsilon(u)+\lambda\,\mathrm{tr}(\epsilon(u))I,\qquad \epsilon(u)=\tfrac12(\nabla u+\nabla u^T),$$3 (Ge et al., 2022).

Initialization is explicit: $$\sigma(u)=2\mu \epsilon(u)+\lambda\,\mathrm{tr}(\epsilon(u))I,\qquad \epsilon(u)=\tfrac12(\nabla u+\nabla u^T),$$4 is the projection of $$\sigma(u)=2\mu \epsilon(u)+\lambda\,\mathrm{tr}(\epsilon(u))I,\qquad \epsilon(u)=\tfrac12(\nabla u+\nabla u^T),$$5, $$\sigma(u)=2\mu \epsilon(u)+\lambda\,\mathrm{tr}(\epsilon(u))I,\qquad \epsilon(u)=\tfrac12(\nabla u+\nabla u^T),$$6 the projection of $$\sigma(u)=2\mu \epsilon(u)+\lambda\,\mathrm{tr}(\epsilon(u))I,\qquad \epsilon(u)=\tfrac12(\nabla u+\nabla u^T),$$7, and $$\sigma(u)=2\mu \epsilon(u)+\lambda\,\mathrm{tr}(\epsilon(u))I,\qquad \epsilon(u)=\tfrac12(\nabla u+\nabla u^T),$$8. Each time step then consists of solving the generalized Stokes problem, updating $$\sigma(u)=2\mu \epsilon(u)+\lambda\,\mathrm{tr}(\epsilon(u))I,\qquad \epsilon(u)=\tfrac12(\nabla u+\nabla u^T),$$9, solving the scalar diffusion problem, and reconstructing $\lambda^\ast,a_0,b_0,\mu,\lambda>0$0 algebraically (Ge et al., 2022).

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 (Wang et al., 22 Mar 2026)	old interface traces in both subproblems	uncoupled variational problems solved in parallel
Robin–Robin splitting (Dalal et al., 2024)	auxiliary Robin data $\lambda^\ast,a_0,b_0,\mu,\lambda>0$1	single decoupled Stokes and Biot solves per time step
Lagrange multiplier partitioned method (Castro et al., 20 Jan 2026)	Schur complement in interface pressure and multipliers	fluid and poroelastic subproblems solved independently at each time step
Multilayered split (Scharf et al., 14 Jul 2025)	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 $\lambda^\ast,a_0,b_0,\mu,\lambda>0$2 and $\lambda^\ast,a_0,b_0,\mu,\lambda>0$3 to enforce tangential-slip and normal-flux/pressure continuity, while all right-hand interface data appear at time level $\lambda^\ast,a_0,b_0,\mu,\lambda>0$4, making the coupling explicit (Wang et al., 22 Mar 2026). 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^\ast,a_0,b_0,\mu,\lambda>0$5, 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^\ast,a_0,b_0,\mu,\lambda>0$6 transfers the latest stress–velocity mismatch to the next time step (Dalal et al., 2024).

In the Lagrange multiplier partitioned method, the coupled weak formulation introduces three interface multipliers,

$\lambda^\ast,a_0,b_0,\mu,\lambda>0$7

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 $\lambda^\ast,a_0,b_0,\mu,\lambda>0$8, after which the fluid velocity, pore pressure, and displacement are recovered by forward solves (Castro et al., 20 Jan 2026).

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 (Scharf et al., 14 Jul 2025).

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 $\lambda^\ast,a_0,b_0,\mu,\lambda>0$9 or $K$0:

$K$1

Here

$K$2

while $K$3 is a cumulative dissipation containing terms such as

$K$4

For $K$5, the method yields a discrete non-increase under the mesh–CFL condition $K$6 (Ge et al., 2022).

For the explicit Stokes–Biot interface split, the analysis introduces a total error energy $K$7 and a dissipative quantity $K$8, then proves the exact energy–dissipation balance

$K$9

After bounding the mixed interface terms $\theta_f>0$0 and the consistency terms $\theta_f>0$1, one obtains a discrete Gronwall inequality and hence unconditional stability in a combined energy-dissipation norm (Wang et al., 22 Mar 2026).

The Robin–Robin scheme has an unconditional stability estimate of the form

$\theta_f>0$2

for zero volume forces, independently of $\theta_f>0$3 and $\theta_f>0$4. The corresponding energy includes fluid kinetic energy, poroelastic kinetic and elastic energy, pore-pressure energy, and an interface contribution involving the auxiliary Robin data $\theta_f>0$5 (Dalal et al., 2024).

By contrast, the multilayered split is proved conditionally stable. With zero external forcing and a CFL-type condition on $\theta_f>0$6 involving mesh size $\theta_f>0$7, thickness $\theta_f>0$8, and physical parameters, the scheme satisfies

$\theta_f>0$9

Under the further parameter constraints $$q:=\nabla\cdot u,\qquad \Phi:=a_0 p+b_0 q,\qquad \Delta:=b_0 p-q-\lambda^\ast q_t,$$0 and $$q:=\nabla\cdot u,\qquad \Phi:=a_0 p+b_0 q,\qquad \Delta:=b_0 p-q-\lambda^\ast q_t,$$1, the paper states a milder CFL condition $$q:=\nabla\cdot u,\qquad \Phi:=a_0 p+b_0 q,\qquad \Delta:=b_0 p-q-\lambda^\ast q_t,$$2 and unconditional stability in the structural coupling parameters (Scharf et al., 14 Jul 2025).

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 $$q:=\nabla\cdot u,\qquad \Phi:=a_0 p+b_0 q,\qquad \Delta:=b_0 p-q-\lambda^\ast q_t,$$3–$$q:=\nabla\cdot u,\qquad \Phi:=a_0 p+b_0 q,\qquad \Delta:=b_0 p-q-\lambda^\ast q_t,$$4–$$q:=\nabla\cdot u,\qquad \Phi:=a_0 p+b_0 q,\qquad \Delta:=b_0 p-q-\lambda^\ast q_t,$$5 elements, the a priori estimate states that, under sufficient regularity, for $$q:=\nabla\cdot u,\qquad \Phi:=a_0 p+b_0 q,\qquad \Delta:=b_0 p-q-\lambda^\ast q_t,$$6 with $$q:=\nabla\cdot u,\qquad \Phi:=a_0 p+b_0 q,\qquad \Delta:=b_0 p-q-\lambda^\ast q_t,$$7 or for $$q:=\nabla\cdot u,\qquad \Phi:=a_0 p+b_0 q,\qquad \Delta:=b_0 p-q-\lambda^\ast q_t,$$8 unconditionally,

$$q:=\nabla\cdot u,\qquad \Phi:=a_0 p+b_0 q,\qquad \Delta:=b_0 p-q-\lambda^\ast q_t,$$9

The same source states that the spatial discretization delivers second-order convergence in $(u,\Delta)$0 for pressures and third-order in $(u,\Delta)$1 for displacement, while backward Euler gives first-order convergence in time (Ge et al., 2022). Numerical tests on $(u,\Delta)$2 report

• $(u,\Delta)$3-error$(u,\Delta)$4 and $(u,\Delta)$5-error$(u,\Delta)$6,
• $(u,\Delta)$7-error$(u,\Delta)$8 and $(u,\Delta)$9-error$p:\Omega\times[0,T]\to\mathbb{R}$00,
• first-order in $p:\Omega\times[0,T]\to\mathbb{R}$01, with representative rates $p:\Omega\times[0,T]\to\mathbb{R}$02, $p:\Omega\times[0,T]\to\mathbb{R}$03, $p:\Omega\times[0,T]\to\mathbb{R}$04, $p:\Omega\times[0,T]\to\mathbb{R}$05, $p:\Omega\times[0,T]\to\mathbb{R}$06, and $p:\Omega\times[0,T]\to\mathbb{R}$07 in the reported manufactured-solution tables (Ge et al., 2022).

For the explicit interface split, the main estimate is

$p:\Omega\times[0,T]\to\mathbb{R}$08

where $p:\Omega\times[0,T]\to\mathbb{R}$09 is the polynomial degree used in the finite element spaces and $p:\Omega\times[0,T]\to\mathbb{R}$10 is the extra regularity exponent in the $p:\Omega\times[0,T]\to\mathbb{R}$11-dual estimates, with $p:\Omega\times[0,T]\to\mathbb{R}$12 on convex domains. The stated interpretation is first-order accuracy in time and optimal spatial convergence rates determined by the approximation spaces (Wang et al., 22 Mar 2026).

For the Lagrange multiplier partitioned method, the convergence theorem gives

$p:\Omega\times[0,T]\to\mathbb{R}$13

with $p:\Omega\times[0,T]\to\mathbb{R}$14 dictated by the polynomial degrees (Castro et al., 20 Jan 2026).

For the Robin–Robin split, the time-discretization error satisfies

$p:\Omega\times[0,T]\to\mathbb{R}$15

and the paper states the time discretization error as $p:\Omega\times[0,T]\to\mathbb{R}$16 (Dalal et al., 2024).

For the multilayered split, the numerical study reports spatial convergence $p:\Omega\times[0,T]\to\mathbb{R}$17 in $p:\Omega\times[0,T]\to\mathbb{R}$18-norm for all variables when using $p:\Omega\times[0,T]\to\mathbb{R}$19 finite elements for poroelastic variables and $p:\Omega\times[0,T]\to\mathbb{R}$20–$p:\Omega\times[0,T]\to\mathbb{R}$21 Taylor–Hood for the fluid, together with temporal convergence $p:\Omega\times[0,T]\to\mathbb{R}$22 for all variables (Scharf et al., 14 Jul 2025). 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 $p:\Omega\times[0,T]\to\mathbb{R}$23 and from which the pore pressure is reconstructed algebraically (Ge et al., 2022). The explicit interface split leads to two uncoupled variational solves per time step, one in $p:\Omega\times[0,T]\to\mathbb{R}$24 and one in $p:\Omega\times[0,T]\to\mathbb{R}$25, with all cross-domain traces lagged to time level $p:\Omega\times[0,T]\to\mathbb{R}$26 (Wang et al., 22 Mar 2026). 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 (Dalal et al., 2024).

The multilayered method is structurally more elaborate. Its Darcy-type substep computes the thick-domain pore pressure $p:\Omega\times[0,T]\to\mathbb{R}$27, plate pressure $p:\Omega\times[0,T]\to\mathbb{R}$28, and Darcy velocities $p:\Omega\times[0,T]\to\mathbb{R}$29 and $p:\Omega\times[0,T]\to\mathbb{R}$30 using fixed-strain data $p:\Omega\times[0,T]\to\mathbb{R}$31, $p:\Omega\times[0,T]\to\mathbb{R}$32, and interface fluid velocity $p:\Omega\times[0,T]\to\mathbb{R}$33. Its structure substep advances $p:\Omega\times[0,T]\to\mathbb{R}$34, $p:\Omega\times[0,T]\to\mathbb{R}$35, $p:\Omega\times[0,T]\to\mathbb{R}$36, $p:\Omega\times[0,T]\to\mathbb{R}$37, and $p:\Omega\times[0,T]\to\mathbb{R}$38 from the updated pressures and interface traction. Its fluid substep then solves the Stokes system for $p:\Omega\times[0,T]\to\mathbb{R}$39 with tangential slip and normal traction prescribed by the new plate–poroelastic data (Scharf et al., 14 Jul 2025). The reported implementation uses FEniCS and standard direct or iterative solvers for each subproblem, together with an additional penalty term

$p:\Omega\times[0,T]\to\mathbb{R}$40

to improve mass-conservation (Scharf et al., 14 Jul 2025).

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:\Omega\times[0,T]\to\mathbb{R}$41. Matrix-vector products with the Schur complement are evaluated through sparse solves with time-invariant matrices $p:\Omega\times[0,T]\to\mathbb{R}$42, $p:\Omega\times[0,T]\to\mathbb{R}$43, and $p:\Omega\times[0,T]\to\mathbb{R}$44, and the preconditioner

$p:\Omega\times[0,T]\to\mathbb{R}$45

is reported to reduce the spectrum of the Schur complement to $p:\Omega\times[0,T]\to\mathbb{R}$46 and to yield iteration counts $p:\Omega\times[0,T]\to\mathbb{R}$47 independent of mesh size (Castro et al., 20 Jan 2026).

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.