Surface Stokes Equations

• Surface Stokes equations are the momentum equations governing viscous, incompressible flow constrained to move tangentially on smooth manifolds.
• They couple surface differential geometry with classical Stokes theory using operators like the Laplace–Beltrami, curvature terms, and symmetric surface strain.
• Recent research employs mixed finite elements, stream-function reductions, and boundary-integral methods to solve these intrinsic PDEs with high accuracy.

Surface Stokes equations are the low-Reynolds-number, surface-incompressible momentum equations for a viscous fluid constrained to move tangentially on a smooth manifold, typically a closed hypersurface embedded in Euclidean space. Their unknowns are a tangential velocity field and a mean-zero surface pressure, and their operator structure couples surface differential geometry with incompressible Stokes theory through the surface gradient, surface divergence, symmetric surface strain, Laplace–Beltrami or Bochner operators, and curvature terms. In current research they appear both as intrinsic geometric PDEs and as limits of bulk thin-film models, and they are treated numerically by mixed finite elements, stream-function reductions, divergence-conforming methods, elliptic reformulations that avoid the discrete inf-sup condition, and boundary-integral formulations (Brandner et al., 2021, Nochetto et al., 18 Aug 2025, Goodwill et al., 23 Feb 2026).

1. Differential-geometric setting

Let $M\subset \mathbb{R}^{d+1}$ be a compact, connected $C^3$ hypersurface without boundary, or in the two-dimensional case $\Gamma\subset \mathbb{R}^3$ a smooth closed surface with unit normal $n$. The orthogonal projector onto the tangent plane is

$P=I-n\otimes n.$

For a scalar $u$ and vector field $v$, the surface differential operators used across the literature include the tangential gradient $\nabla_M u=P\nabla u$, the covariant derivative $\nabla_\Gamma v=P(\nabla v^e)P$, the surface divergence $$\operatorname{div}_\Gamma v=\operatorname{tr}(\nabla_\Gamma v)$$, the scalar Laplace–Beltrami operator $C^3$0, and the Bochner Laplacian on vector fields $C^3$1. The symmetric surface-strain tensor is written either as

$C^3$2

or, equivalently up to notation, $C^3$3 in the finite-element literature (Nochetto et al., 18 Aug 2025, Brandner et al., 2021).

Curvature enters through the Weingarten map $C^3$4 or $C^3$5, the Gaussian curvature $C^3$6, and, in intrinsic formulations, the Ricci tensor. On a two-manifold, $C^3$7. These geometric quantities are not lower-order decoration: they determine equivalences between diffusion operators, affect coercivity, and control mesh-size restrictions in some formulations. The surface-fluid unknown is required to satisfy $C^3$8, so the natural velocity space is a tangential Sobolev space such as

$C^3$9

while the pressure lies in a mean-zero space such as $\Gamma\subset \mathbb{R}^3$0 or $\Gamma\subset \mathbb{R}^3$1 (Nochetto et al., 18 Aug 2025).

2. Governing equations and equivalent operator forms

In primitive variables, one common strong form on a closed surface is

$\Gamma\subset \mathbb{R}^3$2

with tangential force $\Gamma\subset \mathbb{R}^3$3 and tangential velocity $\Gamma\subset \mathbb{R}^3$4 (Brandner et al., 2019). A closely related hypersurface formulation writes

$\Gamma\subset \mathbb{R}^3$5

allowing a source term $\Gamma\subset \mathbb{R}^3$6 and mean-zero pressure (Nochetto et al., 18 Aug 2025). These forms are linked by curvature identities and by the choice of stress operator. In the derivational framework of Brandner–Reusken–Schwering, the Stokes limit $\Gamma\subset \mathbb{R}^3$7 of the tangential surface Navier–Stokes equations yields

$\Gamma\subset \mathbb{R}^3$8

and, using

$\Gamma\subset \mathbb{R}^3$9

one obtains the curvature form

$n$0

for incompressible tangential flow (Brandner et al., 2021).

A rigorous thin-film limit from three-dimensional Navier–Stokes in a curved thin domain leads to the same intrinsic viscous structure. In the special case $n$1 and vanishing boundary friction, Miura obtains

$n$2

which reduces in the Stokes regime to the corresponding steady surface operator with Ricci curvature (Miura, 2020). This connection is important because it identifies the geometric viscous term not merely as a modeling choice, but as a limit of bulk dynamics.

Several numerical works add a zero-order term $n$3 to the momentum equation, for example

$n$4

or

$n$5

This addition is used to avoid global Killing-field issues and regularize the operator on closed surfaces (Hardering et al., 2023, Jankuhn et al., 2020). If $n$6, the derivational literature also records no-slip conditions $n$7 and natural traction conditions $n$8 (Brandner et al., 2021).

3. Variational structure, incompressibility, and Killing fields

The standard weak formulation is a saddle-point problem on a tangential velocity space $n$9 and mean-zero pressure space $P=I-n\otimes n.$0. With

$P=I-n\otimes n.$1

one seeks $P=I-n\otimes n.$2 such that

$P=I-n\otimes n.$3

Equivalent versions replace $P=I-n\otimes n.$4 by curvature-corrected Bochner forms or incorporate source terms in the divergence equation (Krause et al., 2023, Hardering et al., 2023).

Well-posedness relies on a surface Korn inequality and an inf-sup condition. Hardering–Praetorius show coercivity of the velocity bilinear form on $P=I-n\otimes n.$5 and continuous inf-sup stability

$P=I-n\otimes n.$6

while the high-order TraceFEM analysis of Jankuhn–Reusken uses the corresponding surface Korn inequality and inf-sup condition as the continuous foundation for the discrete theory (Hardering et al., 2023, Jankuhn et al., 2020).

A structural complication absent in the flat, simply connected Euclidean setting is the presence of Killing fields. Bonito–Demlow–Licht define a Killing field as a tangential vector field $P=I-n\otimes n.$7 with $P=I-n\otimes n.$8, equivalently an infinitesimal isometry. Their summary records $P=I-n\otimes n.$9 on a generic surface, $u$0 on an axisymmetric surface, and $u$1 on the sphere. Because such modes form a nullspace of the viscous bilinear form, uniqueness requires either orthogonality $u$2, explicit filtering, or the addition of a zero-order regularization term (Bonito et al., 2019).

Tangentiality and incompressibility are logically distinct constraints. Some formulations work directly in tangential spaces, others embed the velocity in $u$3 and enforce $u$4 weakly by a penalty, and divergence-conforming methods enforce $u$5 exactly at the discrete level. This division of roles underlies much of the method design in the recent literature (Bonito et al., 2019, Krause et al., 2023).

4. Reformulations beyond the classical saddle point

On a simply connected oriented surface, every divergence-free tangential field admits a stream function $u$6 such that

$u$7

Substituting this into the Stokes equations yields a fourth-order scalar equation. In the formulation analyzed by Hansbo–Larson–Zahedi and summarized in Brandner et al., the weak form is

$u$8

with

$u$9

Introducing an auxiliary variable $v$0 reduces the problem to a coupled system of two second-order scalar surface PDEs (Brandner et al., 2019, Brandner et al., 2021).

This stream-function route removes the velocity-pressure saddle point entirely, but it requires simple connectivity and then reconstructs $v$1 and $v$2 from scalar solves. For the trace finite element discretization of the coupled $v$3 system, optimal-order bounds are proved, and additional stabilized Laplace–Beltrami problems reconstruct velocity and pressure. On the unit sphere with $v$4, the reported rates are $v$5 for $v$6 and $v$7 in both $v$8- and $v$9-type norms, $\nabla_M u=P\nabla u$0 for $\nabla_M u=P\nabla u$1, $\nabla_M u=P\nabla u$2 for $\nabla_M u=P\nabla u$3, and $\nabla_M u=P\nabla u$4 in the pressure $\nabla_M u=P\nabla u$5-norm with $\nabla_M u=P\nabla u$6 in pressure $\nabla_M u=P\nabla u$7, the last being suboptimal because of geometry error (Brandner et al., 2019).

Nochetto–Shakipov propose a different reformulation on a $\nabla_M u=P\nabla u$8-dimensional $\nabla_M u=P\nabla u$9 hypersurface without boundary. They rewrite the surface Stokes problem as a nonsymmetric indefinite elliptic system governed by two Laplacians, posed for $\nabla_\Gamma v=P(\nabla v^e)P$0. Assuming no geometric error, they prove well-posedness, quasi-best approximation in a robust mesh-dependent $\nabla_\Gamma v=P(\nabla v^e)P$1-norm for any polynomial degree, and optimal $\nabla_\Gamma v=P(\nabla v^e)P$2 error estimates for both velocity and pressure. The key analytical point is that the discrete scheme is stabilized by Gårding inequality, Schatz’s argument, and duality, with a sufficiently small mesh size depending only on the Weingarten map, thereby circumventing the usual discrete inf-sup condition (Nochetto et al., 18 Aug 2025).

This result changes the formulation-dependent narrative around surface Stokes discretization. In the conventional mixed velocity-pressure setting, a Babuška–Brezzi condition is central; in the elliptic reformulation, stability is instead obtained from coercivity up to compact perturbation. The numerical experiments in the same work show equal-order $\nabla_\Gamma v=P(\nabla v^e)P$3 pairs for $\nabla_\Gamma v=P(\nabla v^e)P$4 with observed $\nabla_\Gamma v=P(\nabla v^e)P$5-errors $\nabla_\Gamma v=P(\nabla v^e)P$6 and $\nabla_\Gamma v=P(\nabla v^e)P$7 errors $\nabla_\Gamma v=P(\nabla v^e)P$8, as well as stable mixed-order pairs $\nabla_\Gamma v=P(\nabla v^e)P$9 whose rates are limited by the lower-order field (Nochetto et al., 18 Aug 2025).

5. Finite element discretizations and error theory

A major branch of the literature uses parametric SurfaceFEM on a fitted discrete surface $$\operatorname{div}_\Gamma v=\operatorname{tr}(\nabla_\Gamma v)$$0. Hardering–Praetorius define continuous piecewise-polynomial spaces $$\operatorname{div}_\Gamma v=\operatorname{tr}(\nabla_\Gamma v)$$1 and $$\operatorname{div}_\Gamma v=\operatorname{tr}(\nabla_\Gamma v)$$2, enforce tangentiality by the normal penalty

$$\operatorname{div}_\Gamma v=\operatorname{tr}(\nabla_\Gamma v)$$3

and compare four discrete variants based on two diffusion forms and two divergence forms. They show that the diffusion forms $$\operatorname{div}_\Gamma v=\operatorname{tr}(\nabla_\Gamma v)$$4 and $$\operatorname{div}_\Gamma v=\operatorname{tr}(\nabla_\Gamma v)$$5, and the divergence forms $$\operatorname{div}_\Gamma v=\operatorname{tr}(\nabla_\Gamma v)$$6 and $$\operatorname{div}_\Gamma v=\operatorname{tr}(\nabla_\Gamma v)$$7, are equivalent in the continuous setting but differ discretely by geometric consistency terms such as

$$\operatorname{div}_\Gamma v=\operatorname{tr}(\nabla_\Gamma v)$$8

Using a lifting argument from flat macroelements, they prove discrete inf-sup stability independent of $$\operatorname{div}_\Gamma v=\operatorname{tr}(\nabla_\Gamma v)$$9, and their tangential error estimates yield optimal convergence; on a spherical benchmark with $C^3$00 and Taylor–Hood elements they report fourth-order convergence in the tangential $C^3$01-error (Hardering et al., 2023).

TraceFEM replaces a fitted surface mesh by traces of bulk finite-element spaces on an implicitly or parametrically reconstructed surface. In the higher-order analysis of Jankuhn–Reusken, the discrete method uses Taylor–Hood trace spaces, a penalty term for the normal component, and normal-derivative volume stabilizations

$C^3$02

with $C^3$03, $C^3$04, and $C^3$05. They prove a discrete inf-sup bound uniform in the cut position and an optimal estimate

$C^3$06

with the $C^3$07 term coming solely from geometry and data errors (Jankuhn et al., 2020).

Brandner et al. place parametric TraceFEM and parametric SFEM in a common framework for both velocity-pressure and stream-function formulations. Their reported rates are the expected $C^3$08 in $C^3$09 and $C^3$10 in $C^3$11 for velocity, $C^3$12 in $C^3$13 for pressure, $C^3$14 for tangentiality error $C^3$15, and $C^3$16 for $C^3$17. In their benchmark comparison, SFEM errors are typically $C^3$18–$C^3$19 smaller than TraceFEM for equal degrees of freedom, while both retain optimal slopes (Brandner et al., 2021).

A distinct approach is the divergence-conforming interior-penalty method of Bonito–Demlow–Licht. Using a surface $C^3$20 space and surface Piola mapping, they obtain discrete velocities that are tangential and satisfy $C^3$21 exactly. Interelement $C^3$22-type conformity is imposed weakly by interior-penalty jump terms, and Killing fields are filtered through an auxiliary Stokes eigenproblem. Their error analysis gives

$C^3$23

so that $C^3$24 together with filtering yields optimal $C^3$25-order $C^3$26 (Bonito et al., 2019).

Penalty-only methods remain relevant because of their simplicity. In the fixed-surface framework summarized in the evolving-surface paper of Nestler et al., one replaces the tangential space by $C^3$27 and adds

$C^3$28

with $C^3$29. Under standard ellipticity and inf-sup assumptions, the estimate

$C^3$30

follows, and $C^3$31 yields the expected $C^3$32 rate. The abstract of the same work states that the corresponding evolving-surface Navier–Stokes method exhibits the same optimal order as for stationary surface equations (Krause et al., 2023).

6. Integral equations, derivational links, and methodological themes

An alternative to PDE discretization is the integral-equation formulation of Biliotti–Corona–O’Neil–Rachh. Using two-dimensional Stokeslets in the tangent plane, they represent the solution in terms of a tangential density $C^3$33 and scalar density $C^3$34, obtaining a Fredholm second-kind system

$C^3$35

The operators $C^3$36 are compact on $C^3$37, so the formulation is of “Identity + compact” type and remains well conditioned under refinement. The discretization uses patchwise high-order collocation with Koornwinder polynomials and Vioreanu–Rokhlin nodes, while dense linear algebra is accelerated by proxy-shell compression and recursive skeletonization, leading to $C^3$38–$C^3$39 fast direct solvers (Goodwill et al., 23 Feb 2026).

The numerical behavior reported for this integral method is distinctly high-order. In a compression test, errors below $C^3$40 are achieved with only a few proxy shells. On a slanted torus, fixed collocation order $C^3$41 gives eighth-order convergence in the relative $C^3$42-error for $C^3$43, saturating near $C^3$44, and the direct-solver build time scales essentially linearly in $C^3$45. On an ellipsoid with an added damping term $C^3$46, the relative error is $C^3$47. A star-shaped surface example with prescribed source term $C^3$48 illustrates resolution of complex geometry and local forcing (Goodwill et al., 23 Feb 2026).

The derivational literature clarifies which terms are structural and which are formulation-dependent. Brandner–Reusken–Schwering compare five derivations of evolving-surface Navier–Stokes equations and show that all five yield the same tangential momentum equation. Miura’s thin-film limit then recovers the same viscous operator in a rigorous asymptotic sense. A plausible implication is that the choice between stress-divergence form, Bochner–Ricci form, and curvature-corrected Laplace–Beltrami form is usually a matter of analytical or numerical convenience at the continuous level, provided the corresponding geometric identities are respected (Brandner et al., 2021, Miura, 2020).

Two recurrent misunderstandings are resolved by recent work. First, discrete inf-sup stability is not a universal requirement of the surface Stokes problem itself; it is a requirement of the classical mixed velocity-pressure formulation, whereas the elliptic reformulation of Nochetto–Shakipov obtains stable equal-order and mixed-order discretizations without discrete inf-sup verification when the mesh is sufficiently fine (Nochetto et al., 18 Aug 2025). Second, continuously equivalent diffusion operators are not automatically discretely equivalent in accuracy: Hardering–Praetorius show that the curvature-based form can lose order unless a sufficiently accurate curvature approximation is provided (Hardering et al., 2023).

Taken together, these developments place surface Stokes equations at the intersection of geometric analysis, incompressible flow on manifolds, and high-order scientific computing. The field now contains several analytically mature formulations: classical mixed problems with surface Korn and inf-sup theory, stream-function reductions on simply connected surfaces, divergence-conforming and penalty-based finite elements, inf-sup-free elliptic reformulations, and boundary-integral methods with fast direct solvers. The choice among them is governed less by a single canonical discretization than by geometry representation, treatment of tangentiality and Killing fields, desired conservation properties, and whether the target application favors sparse PDE solvers or dense but high-order integral operators.