Bulk-Surface Stokes System

• The Bulk-Surface Stokes System is a mixed-dimensional framework coupling 3D incompressible bulk flows and 2D surface flows, each with distinct momentum balances and exchange conditions.
• It employs rigorous derivations via variational principles, virtual power, and thermodynamic laws to ensure well-posedness and support advanced numerical schemes.
• Applications span biological membranes, emulsions, and multiphase flows, with numerical methods like TraceFEM, penalty-free FEM, and VEM addressing complex interfacial dynamics.

Searching arXiv for recent and foundational papers on bulk–surface Stokes systems and closely related surface Stokes formulations. Bulk–surface Stokes system denotes a mixed-dimensional viscous flow model in which an incompressible flow in a three-dimensional bulk is coupled to an incompressible viscous flow on an embedded two-dimensional surface or boundary, with the surface treated as a genuine mechanical entity carrying its own momentum balance, tangential kinematics, and, in several formulations, additional phase-field or poroelastic variables. In the Stokes regime, inertia is neglected in both bulk and surface, so the governing subsystem consists of linear momentum balances in the bulk and on the interface, linked by traction, kinematic, and exchange conditions; in diffuse-interface settings this Stokes subsystem is commonly coupled to bulk and surface Cahn–Hilliard equations, while in other mixed-dimensional models it is coupled to surface Biot–Kirchhoff equations or used as the surface component of a broader bulk–surface Navier–Stokes theory (Boschman et al., 2023, Stange, 10 Nov 2025, Knopf et al., 15 Sep 2025).

1. Conceptual scope and geometric setting

In continuum formulations, the bulk is a three-dimensional fluid domain such as $\mathcal{B}\subset\mathbb{R}^3$ or $\Omega\subset\mathbb{R}^3$, and the surface is a two-dimensional manifold such as $\mathcal{S}\subset\mathcal{B}$ or $\Gamma\subset\Omega$. The surface may be closed, may have its own boundary, or may coincide with the boundary of the bulk domain. A recurrent geometric structure is the tangential projector, written either as $$\mathbf{P}_s=\mathbf{I}-\mathbf{n}\otimes\mathbf{n}$$, $T=I_3-n\otimes n$, or $P=I-n\otimes n$, together with surface differential operators $\nabla_s$, $\nabla_\Gamma$, $\operatorname{div}_\Gamma$, and surface symmetric gradients such as $\Omega\subset\mathbb{R}^3$0, $\Omega\subset\mathbb{R}^3$1, or $\Omega\subset\mathbb{R}^3$2 (Boschman et al., 2023, Kone et al., 5 Jun 2026, Goodwill et al., 23 Feb 2026).

A distinctive feature of the bulk–surface viewpoint is that the interface is not merely a geometric boundary. In the finite-thickness continuum theory, the surface is a thin layer with finite thickness $\Omega\subset\mathbb{R}^3$3, homogenized into apparent surface densities and surface stresses living on $\Omega\subset\mathbb{R}^3$4. This introduces the apparent surface density $\Omega\subset\mathbb{R}^3$5, apparent surface free-energy density, and constitutive dependence of surface tension and viscosities on surface state variables. In other formulations, the surface carries its own tangential velocity $\Omega\subset\mathbb{R}^3$6, surface pressure $\Omega\subset\mathbb{R}^3$7, and surface phase field $\Omega\subset\mathbb{R}^3$8, or even plate unknowns such as a normal deflection $\Omega\subset\mathbb{R}^3$9 and a pressure moment $\mathcal{S}\subset\mathcal{B}$0 (Boschman et al., 2023, Stange, 10 Nov 2025, Dassi et al., 4 Aug 2025).

The phrase is not used in a single canonical sense across the literature. Some works study a genuinely coupled 3D–2D fluid system; some analyze a surface Stokes equation that “fits naturally as the surface component of such a coupled bulk–surface system”; and some use a bulk domain only as a computational background, with no feedback from surface flow back into a volumetric flow (Kone et al., 5 Jun 2026, Kaiser et al., 13 Feb 2025).

Literature usage	Defining feature	Representative source
Bulk–surface fluid theory	Bulk and surface each have momentum balance	(Boschman et al., 2023)
Bulk–surface Stokes subsystem	Variable-coefficient bulk and surface Stokes equations on $\mathcal{S}\subset\mathcal{B}$1 and $\mathcal{S}\subset\mathcal{B}$2	(Stange, 10 Nov 2025)
Surface component of bulk–surface flow	Surface Stokes problem posed on $\mathcal{S}\subset\mathcal{B}$3	(Kone et al., 5 Jun 2026)
Mixed 3D–2D coupling beyond surface fluids	Stokes flow in bulk coupled to Biot–Kirchhoff plate on the surface	(Dassi et al., 4 Aug 2025)

2. Governing equations and interface conditions

In the finite-thickness diffuse-interface theory, the bulk Stokes equations are obtained by neglecting inertia in the bulk Navier–Stokes–Cahn–Hilliard-type system. The resulting bulk equations are

$\mathcal{S}\subset\mathcal{B}$4

while the surface Stokes balance is

$\mathcal{S}\subset\mathcal{B}$5

Here $\mathcal{S}\subset\mathcal{B}$6 is a phase-field or Korteweg-type stress, and the surface stress takes a generalized Boussinesq–Scriven form,

$\mathcal{S}\subset\mathcal{B}$7

The coupling is expressed by traction continuity and mixed thermodynamic boundary conditions. In the Stokes regime,

$\mathcal{S}\subset\mathcal{B}$8

and the bulk microflux can satisfy a Robin-type condition

$\mathcal{S}\subset\mathcal{B}$9

In this framework, “bulk–surface Stokes system” therefore means a coupled set of Stokes-type momentum equations in the bulk and on an evolving interface with finite thickness, together with bulk and surface Cahn–Hilliard-type evolution for $\Gamma\subset\Omega$0 and $\Gamma\subset\Omega$1 (Boschman et al., 2023).

In the variable-coefficient two-dimensional boundary-coupled formulation, the stationary bulk–surface Stokes system is written on a bounded $\Gamma\subset\Omega$2 domain $\Gamma\subset\Omega$3 with boundary $\Gamma\subset\Omega$4 as

$\Gamma\subset\Omega$5

$\Gamma\subset\Omega$6

with

$\Gamma\subset\Omega$7

This makes explicit three coupling mechanisms: kinematic coupling through the trace condition, dynamic coupling through the tangential bulk stress $\Gamma\subset\Omega$8, and tangential friction through $\Gamma\subset\Omega$9 (Stange, 10 Nov 2025).

A surface-only Stokes equation is also standard as the surface component of a bulk–surface system. On a smooth closed surface $$\mathbf{P}_s=\mathbf{I}-\mathbf{n}\otimes\mathbf{n}$$0, one form is

$$\mathbf{P}_s=\mathbf{I}-\mathbf{n}\otimes\mathbf{n}$$1

This is the stationary surface Stokes problem analyzed independently in several numerical papers and interpreted there as the natural surface subsystem in a bulk–surface formulation (Kone et al., 5 Jun 2026, Olshanskii et al., 2018).

3. Variational, virtual-power, and thermodynamic derivations

One systematic derivation starts from a principle of virtual powers posed on an arbitrary part $$\mathbf{P}_s=\mathbf{I}-\mathbf{n}\otimes\mathbf{n}$$2, allowing $$\mathbf{P}_s=\mathbf{I}-\mathbf{n}\otimes\mathbf{n}$$3 to be piecewise smooth and permitting edges $$\mathbf{P}_s=\mathbf{I}-\mathbf{n}\otimes\mathbf{n}$$4 where the normal field is discontinuous. The internal virtual power contains a bulk term

$$\mathbf{P}_s=\mathbf{I}-\mathbf{n}\otimes\mathbf{n}$$5

surface terms such as

$$\mathbf{P}_s=\mathbf{I}-\mathbf{n}\otimes\mathbf{n}$$6

and edge contributions like

$$\mathbf{P}_s=\mathbf{I}-\mathbf{n}\otimes\mathbf{n}$$7

After integration by parts and the use of transport theorems on evolving surfaces, one obtains local bulk balances, surface balances, and edge force balance laws. In this setting, jumps of the normal across edges generate jump terms in surface traction, so non-smooth geometry is part of the derivation rather than an afterthought (Boschman et al., 2023).

A second derivational route is based on local mass balance laws, local energy dissipation laws, and the Lagrange multiplier approach. In the bulk–surface Navier–Stokes–Cahn–Hilliard models, this yields bulk and surface chemical potentials, Fick-type fluxes, Newtonian viscous stresses, slip laws, and Robin-type bulk–surface transmission conditions. The total energy combines bulk kinetic energy, bulk Ginzburg–Landau free energy, and surface free energy; the resulting dissipation contains bulk viscous dissipation, surface viscous dissipation, bulk diffusion, surface diffusion, friction, and, when present, chemical exchange terms of the form $$\mathbf{P}_s=\mathbf{I}-\mathbf{n}\otimes\mathbf{n}$$8 (Giorgini et al., 2022, Knopf et al., 15 Sep 2025).

These derivations give the subject a characteristic mixed-dimensional structure. The bulk–surface Stokes subsystem is not introduced ad hoc; it is the linear viscous part selected either from a virtual-power theory with finite surface thickness or from a thermodynamically consistent bulk–surface Navier–Stokes–Cahn–Hilliard system. A plausible implication is that the precise meaning of “surface stress” depends on the parent theory: it may be a Boussinesq–Scriven stress with phase-field corrections, a surface viscous stress acting on a membrane fluid, or the interface term in a Stokes/Biot–Kirchhoff coupling (Boschman et al., 2023, Knopf et al., 15 Sep 2025, Dassi et al., 4 Aug 2025).

4. Analytical structure, well-posedness, and operator theory

For variable coefficients, the bulk–surface Stokes system admits both weak and strong solution theories. In the two-dimensional setting with $$\mathbf{P}_s=\mathbf{I}-\mathbf{n}\otimes\mathbf{n}$$9, $T=I_3-n\otimes n$0, and $T=I_3-n\otimes n$1 uniformly bounded above and below by positive constants, there exists a unique weak solution

$T=I_3-n\otimes n$2

and, for $T=I_3-n\otimes n$3, a unique strong solution

$T=I_3-n\otimes n$4

with

$T=I_3-n\otimes n$5

The proof combines a bulk-surface Korn inequality, elliptic regularity for the bulk and surface subproblems, pressure reconstruction, and a Schauder fixed-point argument for the traction-coupled variable-coefficient system (Stange, 10 Nov 2025).

In the constant-coefficient setting, a bulk-surface Stokes operator $T=I_3-n\otimes n$6 can be defined on the divergence-free product space. It is positive, self-adjoint on $T=I_3-n\otimes n$7 with compact inverse, and therefore admits a countable sequence of positive eigenvalues with eigenfunctions forming an orthonormal basis of $T=I_3-n\otimes n$8. Those eigenfunctions are used as the Galerkin basis in the existence proof for global weak solutions to a thermodynamically consistent bulk–surface Navier–Stokes–Cahn–Hilliard model with an additional surface Navier–Stokes equation (Knopf et al., 15 Sep 2025).

The analytical role of the bulk–surface Stokes system is therefore twofold. First, it is the linear elliptic subsystem needed for higher-order regularity and uniqueness arguments in bulk–surface Navier–Stokes–Cahn–Hilliard models. Second, it is an operator-theoretic backbone for semi-Galerkin approximations, spectral decompositions, and compactness arguments. This also explains why some papers treat “surface Stokes” separately: the surface operator is often the difficult geometric part of the coupled mixed-dimensional system (Stange, 10 Nov 2025, Knopf et al., 15 Sep 2025).

5. Numerical formulations and discretisation strategies

The numerical literature is heterogeneous because the surface constraint $T=I_3-n\otimes n$9, incompressibility, curvature, and mixed-dimensional coupling are difficult to enforce simultaneously. Trace finite element methods approximate the surface problem on a background bulk mesh. A low-order TraceFEM for the surface Stokes problem uses $P=I-n\otimes n$0 bulk finite elements for both the velocity and the pressure, introduces a penalty term to enforce tangentiality, uses additional normal-derivative volume stabilization for velocity and pressure, and proves stability and optimal order error bounds in the surface $P=I-n\otimes n$1 and $P=I-n\otimes n$2 norms, with an $P=I-n\otimes n$3 geometric consistency error (Olshanskii et al., 2018). Higher-order parametric trace FEMs employ generalized Taylor–Hood pairs on tetrahedral bulk meshes, quantify geometric errors from approximate parametric surface representation, and prove optimal order convergence in the isoparametric regime (Jankuhn et al., 2020, Jankuhn et al., 2019).

A different line of work removes the penalty. A tangential and penalty-free finite element method constructs surface MINI spaces whose velocity fields are tangential, are not $P=I-n\otimes n$4-conforming, but do lie in $P=I-n\otimes n$5, and do not require penalization to achieve optimal convergence rates; numerical experiments also indicate optimal-order convergence of nonconforming tangential surface Taylor–Hood $P=I-n\otimes n$6 elements (Demlow et al., 2023). Another development constructs an exactly divergence-free Scott–Vogelius finite element method for the surface Stokes problem on curved Clough–Tocher triangulations; the method simultaneously enforces tangentiality and incompressibility exactly, proves inf-sup stability, and derives optimal-order convergence in the isoparametric regime (Kone et al., 5 Jun 2026). A Lagrange-multiplier surface finite element method based on Taylor–Hood elements derives a new inf-sup condition for the additional multiplier and proves optimal velocity convergence in energy and tangential $P=I-n\otimes n$7 norms, along with optimal $P=I-n\otimes n$8 convergence for the two pressures (Elliott et al., 2024).

Integral-equation methods provide a distinct route. A parametrix for the surface Stokes equation represents the surface operator as a Fredholm integral equation of the second kind built from two-dimensional Stokeslets, with tangentiality built into the ansatz and a proxy shell method for fast direct solvers of the resulting dense linear systems (Goodwill et al., 23 Feb 2026).

For genuinely coupled mixed-dimensional systems, virtual element methods have also been developed. In the Stokes/Biot–Kirchhoff bulk–surface model, a stable VEM is established for the monolithic 3D–2D problem, with a discrete inf-sup condition under a small mesh assumption through a Fortin interpolant that requires only $P=I-n\otimes n$9-regularity for the Stokes problem; the well-posedness of the monolithic discrete formulation and the optimal convergence in the energy norm are proved theoretically and confirmed numerically (Dassi et al., 4 Aug 2025).

Discretisation family	Characteristic feature	Representative source
TraceFEM with penalty	Background bulk mesh, tangentiality via penalty	(Olshanskii et al., 2018)
Higher-order trace FEM	Parametric geometry, Taylor–Hood pairs	(Jankuhn et al., 2020)
Penalty-free tangential FEM	Tangential velocity spaces in $\nabla_s$0	(Demlow et al., 2023)
Exactly divergence-free surface FE	Scott–Vogelius on curved Clough–Tocher meshes	(Kone et al., 5 Jun 2026)
Lagrange-multiplier SFEM	Tangentiality enforced by multiplier	(Elliott et al., 2024)
Integral equations	Fredholm second-kind formulation, fast direct solver	(Goodwill et al., 23 Feb 2026)
Monolithic VEM for 3D–2D coupling	Bulk Stokes + surface Biot–Kirchhoff	(Dassi et al., 4 Aug 2025)

6. Applications, related models, and interpretive issues

The motivating applications are interfacial and membrane flows with internal structure. The finite-thickness continuum theory is motivated by biological membranes and cell cortex flows, two-phase flows with interfacial regions such as emulsions and polymer blends, and thin layers with phase segregation on substrates. In these settings, the finite surface thickness and apparent surface densities allow in-plane diffusion and phase separation on the interface, while bulk–surface coupling handles mass exchange and interfacial rheology in a thermodynamically consistent way (Boschman et al., 2023).

Biological motivation is particularly explicit in the thermodynamically consistent viscous fluid-mixture models with a surface Navier–Stokes equation. There the inclusion of surface hydrodynamics is linked to the fluid mosaic model of cell membranes, in which the surface of biological cells is interpreted as a thin layer of viscous fluids. This suggests a surface equation with its own surface viscosity, surface inertia, and surface phase segregation rather than a purely passive dynamic boundary condition (Knopf et al., 15 Sep 2025).

The phrase also appears in broader mixed-dimensional settings. A coupled 3D–2D model with a free fluid governed by Stokes flow in the bulk and a poroelastic plate described by the Biot–Kirchhoff equations on the surface is analyzed and discretized by VEM, with an application to immune isolation using encapsulation with silicon nanopore membranes (Dassi et al., 4 Aug 2025). This suggests that, in the broader PDE literature, “bulk–surface Stokes system” can denote a Stokes-governed bulk coupled not only to a surface fluid but to other lower-dimensional physics on the surface.

A common misconception is that the phrase always denotes a single canonical PDE system. The literature instead contains at least three stable meanings. In one, it is the Stokes limit of a bulk–surface Navier–Stokes–Cahn–Hilliard model with finite surface thickness and phase segregation. In another, it is a variable-coefficient coupled bulk-and-boundary Stokes problem with trace coupling, tangential stress, and friction. In a third, especially in several numerical papers, the bulk may be only a computational background mesh or a co-area integration domain rather than a physical fluid region (Boschman et al., 2023, Stange, 10 Nov 2025, Kaiser et al., 13 Feb 2025).

That distinction matters computationally and conceptually. In “Bulk Trace FEM,” for example, the bulk domain is a computational background and there is no feedback from surface flow back into a volumetric flow; by contrast, in the coupled viscous-mixture and Stokes/Biot–Kirchhoff models, bulk and surface exchange momentum through traction or stress terms and are part of a single monolithic mixed-dimensional problem (Kaiser et al., 13 Feb 2025, Dassi et al., 4 Aug 2025). In electrokinetic transport, an analogous bulk–surface structure appears when the bulk Stokes flow in a microchannel is coupled to diffuse electric-double-layer surface currents and effective slip conditions at charged walls; the reported result that, without a surface current, bulk advection is strongly suppressed illustrates how a surface transport mechanism can control bulk Stokes-scale motion (Nielsen et al., 2014).

The modern subject is therefore best understood as a family of rigorously formulated mixed-dimensional viscous systems. Their common core is the coupling of a bulk incompressible Stokes-type balance to a lower-dimensional surface balance through geometry, traction, kinematics, and, in many models, thermodynamically derived exchange laws.