1. Introduction
When surface dissipation dominates bulk dissipation (i.e. when the Saffman–Delbrück length far exceeds the system size), the viscous flow of evolving fluid films can be modelled without accounting for the bulk fluid (Saffman & Delbrück Reference Saffman and Delbrück1975; Arroyo & DeSimone Reference Arroyo and DeSimone2009). The governing equations, known as the evolving Stokes equations, have not only been of theoretical interest, but also found practical applications in engineering, such as the study of foam (Exerowa & Kruglyakov Reference Exerowa and Kruglyakov1997). Their use also has a long-standing history in biophysical contexts, addressing fundamental problems across scales – from modelling subcellular structures such as lipid bilayers to tissue-level phenomena such as epithelial monolayers (Al-Izzi & Morris Reference Al-Izzi and Morris2021).
Historically, the governing equations were formulated by Scriven (Reference Scriven1960) to study foam instability. When coupled with bending elasticity, the evolving Stokes equations serve as a common model for lipid membranes (Arroyo & DeSimone Reference Arroyo and DeSimone2009). Recently, there has been substantial interest in modelling cell and tissue growth by coupling the viscous layer with additional anisotropy and active processes, such as those seen in the cellular cortex and active nematic fluids (Metselaar, Yeomans & Doostmohammadi Reference Metselaar, Yeomans and Doostmohammadi2019; Torres-Sánchez, Millán & Arroyo Reference Torres-Sánchez, Millán and Arroyo2019; Al-Izzi & Morris Reference Al-Izzi and Morris2023). The evolving Stokes equations, along with extensions to the full Navier–Stokes equations, have been derived independently through various principles (Scriven Reference Scriven1960; Arroyo & DeSimone Reference Arroyo and DeSimone2009; Koba, Liu & Giga Reference Koba, Liu and Giga2017; Miura Reference Miura2018; Torres-Sánchez et al. Reference Torres-Sánchez, Millán and Arroyo2019; Reuther, Nitschke & Voigt Reference Reuther, Nitschke and Voigt2020; Sahu et al. Reference Sahu, Omar, Sauer and Mandadapu2020; Al-Izzi & Morris Reference Al-Izzi and Morris2023) and are compared in Brandner, Reusken & Schwering (Reference Brandner, Reusken and Schwering2022). Our approach is closest to the formulation of Arroyo & DeSimone (Reference Arroyo and DeSimone2009) via Onsager's variational principle for dissipative systems.
To computationally solve the evolving Stokes equations, most existing methods tackle the partial differential equations by explicitly splitting them into normal and tangential components. This involves a vector-valued equation for the tangent velocity and a scalar equation for the normal velocity on a manifold. When viewed this way, the system of equations is challenging to solve. First, the tangent equation is tensor-valued on an evolving Riemannian manifold, necessitating specialized techniques for covariant differentiation (Knöppel et al. Reference Knöppel, Crane, Pinkall and Schröder2013; Gross & Atzberger Reference Gross and Atzberger2018; Nestler, Nitschke & Voigt Reference Nestler, Nitschke and Voigt2019; Voigt Reference Voigt2019; Torres-Sánchez, Santos-Oliván & Arroyo Reference Torres-Sánchez, Santos-Oliván and Arroyo2020). Some methods avoid the tensor-valued equation of tangent velocity through streamfunction formulations, although this approach is limited to simply connected domains (Torres-Sánchez et al. Reference Torres-Sánchez, Millán and Arroyo2019) unless additional topological techniques are employed (Yin et al. Reference Yin, Nabizadeh, Wu, Wang and Chern2023). Second, the equations in the decomposed form are explicitly coupled with surface curvatures, which, in classic numerical approximations, require representations using high-order basis functions (Sahu et al. Reference Sahu, Omar, Sauer and Mandadapu2020; Krause & Voigt Reference Krause and Voigt2023).
In this study, instead of directly tackling the classic equations, we return to the fundamental governing kinematics and principles of a viscous surface film. We replace the system of equations expressed in tangent–normal splitting with a simple, coordinate-free differential-geometric formulation. We show that our resulting equations agree with earlier formulations by Scriven (Reference Scriven1960) and Arroyo & DeSimone (Reference Arroyo and DeSimone2009) but in a more elegant form. Our approach also directly carries over to the discrete setting. The abstraction directly leads to a discretization for the strain rate tensor on discrete meshes. We then construct a discrete Rayleighian for the system and derive the discrete Stokes equations via Onsager's variational principle, leading to a simple linear system. We provide a self-contained exposition of both the continuous and discrete models, as well as the numerical methods used to solve the system of equations on arbitrary geometries and topologies.
2. Theory
We derive the equations of motion for an evolving surface driven by viscous Stokes flow using Onsager's variational principle. Specifically, a dissipation functional, called the Rayleighian, is assigned to each position and velocity state of a deformable surface. The equations of motion are then obtained by finding the velocity that minimizes this functional.
We adopt a differential-geometric approach, similar to Marsden & Hughes (Reference Marsden and Hughes1994), to describe continuum mechanics. Tensors are expressed using fibre bundle notation for clarity. Readers unfamiliar with this terminology can refer to § 2.1 for a brief overview. We first define relevant tensors and differential operators in § 2.2, followed by the definition of the Rayleighian for the system in § 2.3. For simplicity in presenting our smooth theory, all functions and tensors are assumed to be of 
$C^\infty$ class.
2.1. Terminology and notation
On a manifold 
$M$, the ‘tangent space’ 
$T_pM$ is a vector space that provides a local linear approximation of the manifold at 
$p \in M$. The disjoint union of all tangent spaces forms the ‘tangent bundle’ 
$TM = \bigcup _p T_pM$, where each 
$T_pM$ is also called the ‘fibre’ of the tangent bundle 
$TM$ at the ‘base point’ 
$p\in M$ on its ‘base manifold’ 
$M$. The dual bundle of 
$TM$ is the ‘cotangent bundle’, 
$T^*M$, whose fibres are the dual spaces of the tangent spaces. The 
$(m, n)$-typed tensor bundle is written as 
$TM^{\otimes m, n} = TM^{\otimes m} \otimes T^*M^{\otimes n}$. A tensor field 
$\boldsymbol {\varPsi }$ on 
$M$ is a ‘section’ of the tensor bundle, denoted as 
$\boldsymbol {\varPsi }\in \varGamma (TM^{\otimes m, n})$, which assigns a tensor 
$\boldsymbol {\varPsi }|_p\in T_pM^{\otimes m, n}$ to each point 
$p\in M$. When a basis section is chosen, its components, 
$\varPsi _{i_1 \ldots i_n}^{j_1 \ldots j_m}$, have 
$m$ contravariant indices and 
$n$ covariant indices.
For a map 
$\varphi : M \rightarrow W$, where 
$W$ is another manifold, the ‘pullback bundle’ 
$T_{\varphi }W = \varphi ^*TW$ is defined with 
$M$ as its base manifold, while the fibre at 
$p\in M$, 
$(T_\varphi W)|_p = T_{\varphi (p)}W$, is the tangent space to 
$W$ at 
$\varphi (p) \in W$.
Note that tensors can be identified with linear maps; for example, 
$\boldsymbol {\varPsi }|_p \in T_pM^{\otimes m, n}$ acts as a linear map 
$\boldsymbol {\varPsi }|_p: T_pM \rightarrow T_pM^{\otimes m, n-1}$.
2.2. Differential geometry of an evolving surface
We describe a deformable surface as follows. Let 
$M$ be a 2-dimensional closed manifold (
$\partial M = \varnothing$) with arbitrary genus, representing the ‘material/Lagrangian space’. Let 
$W$ be a three-dimensional (3-D) Riemannian manifold representing the ‘world/Eulerian space’. Let 
${{\boldsymbol{\mathsf{g}}}}\in \varGamma (T^*W\otimes _{symm}T^*W)$ (i.e. a symmetric (0, 2) tensor field) denote the Riemannian metric tensor on 
$W$. The common setup is to set 
$W = \mathbb{R}^3$, the Euclidean 3-space, with 
$\boldsymbol{\mathsf{g}}[\![ \boldsymbol {V},\boldsymbol {W} ]\!] = \boldsymbol {V} \boldsymbol {\cdot } \boldsymbol {W}$ being the Euclidean inner product for any 
$\boldsymbol {V},\boldsymbol {W}\in \mathbb{R}^3$. However, our theory is not limited to this assumption. The ‘position’ of a deformable surface is an embedding function 
$\varphi \colon M\hookrightarrow W$ from the material space to the world space. A time-dependent deformable surface, or an evolving surface, is a time-dependent embedding 
$\varphi (t) = \varphi _t\colon M\hookrightarrow W$. The ‘material velocity’ of an evolving surface 
$\varphi (t)$ is given by 
$\dot \varphi _t := \partial _t\varphi _t$, which is a ‘pullback’ vector field 
$\dot \varphi _t\in \varGamma (T_{\varphi _t}W)$ with its base point in the material frame 
$M$ (cf. § 2.1). (The velocity 
$\dot \varphi _t|_p$ is evaluated at each point 
$p\in M$ in the material domain, and its value 
$\dot \varphi _t|_p\in T_{\varphi _t(p)}W$ is a tangent vector to the world space based at 
$\varphi _t(p)\in W$.)
Each material velocity 
$\dot \varphi _t\in \varGamma (T_{\varphi _t}W)$ can be extended to a ‘spatial velocity’ over 
$W$, denoted by 
$\boldsymbol {U}_t\in \varGamma (TW)$, such that 
$\dot \varphi _t = \boldsymbol {U}_t\circ \varphi _t$, or equivalently, 
$\dot \varphi _t|_p = \boldsymbol {U}_t|_{\varphi _t (p)}$ for 
$\varphi _t(p) \in W$. Note that this extension is not unique, as the assignment of 
$\boldsymbol {U}_t$ at locations away from the surface 
$\varphi _t(M)\subset W$ is arbitrary.
2.2.1. Deformation gradient and Cauchy–Green tensor
The differential of the positioning 
$\varphi \colon M\hookrightarrow W$ of a surface is denoted by 
${{\boldsymbol{\mathsf{F}}}}:= d\varphi \in \varGamma (T_{\varphi }W\otimes T^*M)$. This differential is also known as the pushforward map, the tangent map, or the ‘deformation gradient’. At each point 
$p\in M$, the two-point tensor 
${{\boldsymbol{\mathsf{F}}}}|_p$ can be identified with a linear map (cf. § 2.1) 
${{\boldsymbol{\mathsf{F}}}}|_p = d\varphi |_p\colon T_pM\rightarrow T_{\varphi (p)}W$ that realizes a material tangent vector as a 3-D world vector, as expected for the deformation gradient. (At each point 
$p\in M$, the tensor 
${{\boldsymbol{\mathsf{F}}}}|_p\in T_{\varphi (p)}W\otimes T_{p}^*M$ pairs with a vector in 
$T_pM$ and returns a value in 
$T_{\varphi (p)}W$. Since it relates quantities in two different spaces 
$M$ and 
$W$, it is called a two-point tensor.)
The embedding 
$\varphi \colon M\hookrightarrow W$ induces a Riemannian metric (the first fundamental form) 
$\boldsymbol{\mathsf{I}}\in \varGamma (T^*M\otimes _{symm}T^*M)$ on 
$M$ from the Riemannian structure 
${{\boldsymbol{\mathsf{g}}}}$ on 
$W$. This induced metric is defined by 
$\boldsymbol{\mathsf{I}}[\![\boldsymbol {v},\boldsymbol {w}]\!]:={{\boldsymbol{\mathsf{g}}}}[\![ {{\boldsymbol{\mathsf{F}}}}\boldsymbol {v},{{\boldsymbol{\mathsf{F}}}}\boldsymbol {w}]\!] = \boldsymbol {v}^\top {{\boldsymbol{\mathsf{F}}}}^\top {{\boldsymbol{\mathsf{g}}}}{{\boldsymbol{\mathsf{F}}}}\boldsymbol {w}$ for each 
$\boldsymbol {v},\boldsymbol {w}\in T_pM$, 
$p\in M$. In continuum mechanics, 
$\boldsymbol{\mathsf{I}} = {{\boldsymbol{\mathsf{F}}}}^\top {{\boldsymbol{\mathsf{g}}}}{{\boldsymbol{\mathsf{F}}}}$ is known as the ‘right Cauchy–Green tensor’. Here we identify tensors 
$\boldsymbol{\mathsf{I}}, \boldsymbol{\mathsf{g}}, \boldsymbol{\mathsf{F}}$ as linear maps 
$\boldsymbol{\mathsf{I}}_p \colon T_pM \rightarrow T_p^*M$, 
$\boldsymbol{\mathsf{g}}_q \colon T_qW \rightarrow T_q^*W$, 
$\boldsymbol{\mathsf{F}}_p\colon T_pM\rightarrow T_{\varphi(p)}W$. Then we have the compositional relation 
$\boldsymbol{\mathsf{I}} = {{\boldsymbol{\mathsf{F}}}}^\top {{\boldsymbol{\mathsf{g}}}}{{\boldsymbol{\mathsf{F}}}}$, where 
${{\boldsymbol{\mathsf{F}}}}_p^\top \colon T_{\varphi (p)}^*W\rightarrow T_p^*M$ is the adjoint of 
${{\boldsymbol{\mathsf{F}}}}_p$.)
For brevity, we denote the induced metric as 
$\langle {\cdot },{\cdot }\rangle = \boldsymbol{\mathsf{I}}[\![{\cdot },{\cdot }]\!]$ and the 
$L_2$ inner product as 
$\langle \!\langle {\cdot },{\cdot }\rangle \!\rangle = \int _M\langle {\cdot },{\cdot }\rangle \, \mathrm {d}A$ where 
$\mathrm {d} A$ is the area form from the metric 
$\boldsymbol{\mathsf{I}}$. These inner products define the vector norm 
$|{\cdot }|^2$ and the 
$L_2$ norm 
$\|{\cdot }\|^2$. Dual pairings between primal (contravariant) and dual (covariant) tensors are denoted by 
$\langle {\cdot } \,|\, {\cdot } \rangle$ and 
$\langle \!\langle {\cdot } \,|\, {\cdot } \rangle \!\rangle = \int _M \langle {\cdot } \,|\, {\cdot } \rangle \,\mathrm {d} A$. The 
$L_2$ inner product between two scalar functions 
$f$ and 
$g$ is denoted by 
$\langle \!\langle\, f, g\rangle \!\rangle = \int _{M}fg\, \mathrm {d} A$.
2.2.2. Strain rate tensor
For an evolving surface 
$\varphi _t$, the induced metric 
$\boldsymbol{\mathsf{I}}_t\in \varGamma (T^*M\otimes _{symm}T^*M)$ varies with time 
$t$. The ‘strain rate tensor’ 
${{\boldsymbol{\mathsf{E}}}}_t\in \varGamma (T^*M\otimes _{symm}T^*M)$ is defined as half the rate of change of the induced metric: 
$2{{\boldsymbol{\mathsf{E}}}}_t := \dot{\boldsymbol{\mathsf{I}}}_t$.
To compute the time derivative of 
$\boldsymbol{\mathsf{I}}_t = {{\boldsymbol{\mathsf{F}}}}^\top _t{{\boldsymbol{\mathsf{g}}}}{{\boldsymbol{\mathsf{F}}}}_t$ we first calculate the time derivative of 
${{\boldsymbol{\mathsf{F}}}}_t$. Let 
$\boldsymbol {U}_t\in \varGamma (TW)$ be any extended velocity field such that 
$\dot \varphi _t = \boldsymbol {U}_t\circ \varphi _t$. Then,
\begin{equation} \dot{{{\boldsymbol{\mathsf{F}}}}} = d \dot \varphi_t = d (\boldsymbol{U}_t \circ \varphi_t) = {\boldsymbol{\nabla}} (\boldsymbol{U}_t \circ \varphi_t) \circ d \varphi_t \equiv (\boldsymbol{\nabla} \boldsymbol{U}_t) {{\boldsymbol{\mathsf{F}}}}, \end{equation}
where 
${\boldsymbol {\nabla }: \varGamma (TW) \rightarrow \varGamma (TW\otimes T^*W)}$ is the Levi-Civita covariant derivative compatible with the metric 
${{\boldsymbol{\mathsf{g}}}}$. When applied to a vector 
$\boldsymbol {v}\in \varGamma (TM)$, 
$\dot{{{\boldsymbol{\mathsf{F}}}}}_t\boldsymbol {v} = {\boldsymbol {\nabla }}_{{{\boldsymbol{\mathsf{F}}}}_t\boldsymbol {v}}\boldsymbol {U}_t$, which is a directional derivative of 
$\boldsymbol {U}_t$ tangential to the surface. In particular, (2.1) depends only on the values of 
$\boldsymbol {U}_t\circ \varphi _t$ at the surface and not on the choice of its extension 
$\boldsymbol {U}_t$ over 
$W$.
With a time-invariant metric 
${{\boldsymbol{\mathsf{g}}}}$ and (2.1), the strain rate tensor is given by
\begin{equation} 2{{\boldsymbol{\mathsf{E}}}}_t = \dot{\boldsymbol{\mathsf{I}}}_t = \dot{{{\boldsymbol{\mathsf{F}}}}}_t^\top {{\boldsymbol{\mathsf{g}}}} {{\boldsymbol{\mathsf{F}}}}_t + {{\boldsymbol{\mathsf{F}}}}_t^\top {{\boldsymbol{\mathsf{g}}}} \dot{{{\boldsymbol{\mathsf{F}}}}}_t = {{\boldsymbol{\mathsf{F}}}}_t^\top (({\boldsymbol{\nabla}} \boldsymbol{U}_t)^\top {{\boldsymbol{\mathsf{g}}}} + {{\boldsymbol{\mathsf{g}}}} {\boldsymbol{\nabla}} \boldsymbol{U}_t ) {{\boldsymbol{\mathsf{F}}}}_t. \end{equation}
When applied to vectors, 
$2{{\boldsymbol{\mathsf{E}}}}_t [\![ \boldsymbol {v}, \boldsymbol {w} ]\!]|_p = ({\boldsymbol {\nabla }}_{\boldsymbol {F}_t \boldsymbol {v}} \boldsymbol {U}_t)\boldsymbol {\cdot }(\boldsymbol {F}_t \boldsymbol {w}) + ({\boldsymbol {\nabla }}_{\boldsymbol {F}_t \boldsymbol {w}} \boldsymbol {U}_t)\boldsymbol {\cdot }(\boldsymbol {F}_t \boldsymbol {v})$ for 
$\boldsymbol {v}, \boldsymbol {w} \in \varGamma (T_pM)$. Again, 
${{\boldsymbol{\mathsf{E}}}}_t$ depends only on the values of 
$\boldsymbol {U}_t\circ \varphi _t$ and not on the choice of its extension 
$\boldsymbol {U}_t$.
2.2.3. Differential operators
The operator that takes each velocity field 
$\dot \varphi = \boldsymbol {U}\circ \varphi$ at the surface to its induced strain rate tensor (2.2) is called the Killing operator:
\begin{equation} {\mathcal{K}}: \varGamma(T_{\varphi}W) \rightarrow \varGamma(T^*M \otimes_{symm} T^*M), \quad {\mathcal{K}} \dot \varphi := {{\boldsymbol{\mathsf{E}}}} = {{\boldsymbol{\mathsf{F}}}}^\top \frac{({\boldsymbol{\nabla}}\boldsymbol{U})^\top{{\boldsymbol{\mathsf{g}}}} + {{\boldsymbol{\mathsf{g}}}}{\boldsymbol{\nabla}}\boldsymbol{U}}{2}{{\boldsymbol{\mathsf{F}}}}. \end{equation}The trace of the strain rate tensor defines the ‘divergence operator’
\begin{equation} \operatorname{DIV} \colon \varGamma(T_\varphi W) \rightarrow C^{\infty}(M),\quad \operatorname{DIV} \dot \varphi := \operatorname{tr} ({\mathcal{K}} \dot \varphi), \end{equation}which measures the rate of change of area induced by the given velocity field. Define the ‘gradient operator’ as the negative adjoint of the divergence operator:
\begin{equation} \operatorname{GRAD} ={-} \operatorname{DIV}^*: C^{\infty}(M) \rightarrow \varGamma(T_\varphi^*W),\quad \langle\!\langle\operatorname{GRAD} f\,|\,\dot \varphi\rangle\!\rangle ={-} \langle\!\langle\, f, \operatorname{DIV}\dot \varphi\rangle\!\rangle, \end{equation}
for all 
$f\in C^\infty (M)$. The ‘covariant divergence’ is the negative adjoint of the Killing operator:
\begin{align} \operatorname{DIV}^{\boldsymbol{\nabla}} ={-} {\mathcal{K}}^*: \varGamma(TM \otimes_{\operatorname{symm}} TM) \rightarrow \varGamma(T_\varphi^*W),\quad \langle\!\langle \operatorname{DIV}^{\boldsymbol{\nabla}} {{\boldsymbol{\mathsf{T}}}}\,|\,\dot \varphi\rangle\!\rangle ={-}\langle\!\langle {{\boldsymbol{\mathsf{T}}}}\,|\,{\mathcal{K}}\dot \varphi\rangle\!\rangle, \end{align}
for all symmetric tensor fields 
${{\boldsymbol{\mathsf{T}}}}\in \varGamma (TM\otimes _{symm}TM)$. Here the dual pairing 
$\langle {{\boldsymbol{\mathsf{T}}}}\,|\,{{\boldsymbol{\mathsf{E}}}}\rangle$ between 
$(0,2)$ and 
$(2,0)$ tensors is 
$\langle {{\boldsymbol{\mathsf{T}}}}\,|\,{{\boldsymbol{\mathsf{E}}}}\rangle = \operatorname {tr}({{\boldsymbol{\mathsf{T}}}}^\top {{\boldsymbol{\mathsf{E}}}})= \sum _{ij} T^{ij}E_{ij}$.
Notably, the operators 
${\mathcal {K}}$, 
$\operatorname {DIV}$, 
$\operatorname {GRAD}$ and 
$\operatorname {DIV}^{\boldsymbol {\nabla }}$, derived from the strain rate tensor 
${{\boldsymbol{\mathsf{E}}}}$ of the evolving surface, differ from the commonly used intrinsic operators 
$\mathscr {K}$, 
$\operatorname {div}$, 
$\operatorname {grad}$ and 
$\operatorname {div}^{\boldsymbol {\nabla }}$. As shown through the normal–tangent decomposition in the Appendix, these differences arise from additional contributions related to curvature (cf. A3–A5).
2.3. Onsager's variational principle for an evolving viscous fluid film
The relaxation dynamics of an evolving viscous fluid film is governed by Onsager's variational principle, which states that Stokes flow follows the dynamics of the least energy dissipation (Arroyo & DeSimone Reference Arroyo and DeSimone2009; Doi Reference Doi2011). The dissipation functional, or Rayleighian, captures the system's dissipation rate, and the governing flow is the stationary condition of this functional subject to the incompressibility constraint.
Consider a viscous dynamical system where the dissipation functional 
${\mathcal {D}}$ comprises the rate of viscous dissipation, and power inputs from an external force 
$\boldsymbol {B}\in \varGamma (T^*W)$ and a shape-dependent potential energy 
$V_\varphi$:
\begin{align} {\mathcal{D}}(\boldsymbol{U}) &= \frac{1}{2} \langle\!\langle {{\boldsymbol{\mathsf{E}}}} \,|\, {{{\boldsymbol{\mathsf{T}}}}} \rangle\!\rangle -\langle\!\langle \boldsymbol{B} \,|\, \boldsymbol{U} \rangle\!\rangle + \dot V_\varphi\nonumber\\ & = \frac{1}{2} \int_M \langle {{\boldsymbol{\mathsf{E}}}} \,|\, {{\boldsymbol{\mathsf{T}}}} \rangle \,\mathrm{d} A - \int_{\varphi(M)} \langle \boldsymbol{B} \,|\, \boldsymbol{U} \rangle \,\mathrm{d} A + \int_M \langle \delta_\varphi V_\varphi\,|\, \dot \varphi \rangle \,\mathrm{d}A. \end{align}
Here, 
$\boldsymbol{\mathsf{T}} = 2 \mu \boldsymbol{\mathsf{E}} \in \varGamma(TM \otimes_{{\textit{symm}}} TM)$ is the viscous stress tensor (assumed Newtonian with constant viscosity 
$\mu$), 
$\delta _\varphi V_\varphi$ is the conservative force obtained by taking the variation of 
$V_\varphi$, and 
$\dot \varphi = \boldsymbol {U} \circ \varphi$ is the material velocity.
Onsager's variational principle determines the Stokes flow by solving the constrained optimization problem:
\begin{equation} \min_{\boldsymbol{U}}{\mathcal{D}}(\boldsymbol{U}) \quad \text{subject to} \ \operatorname{DIV} \boldsymbol{U} = 0. \end{equation}
Using the fluid pressure 
$p \in C^\infty (M)$ as a Lagrange multiplier, this problem can be formulated as a minimax problem involving the augmented Rayleighian 
${\mathcal {R}}$:
\begin{equation} \min_{\boldsymbol{U}} \max_{p} {\mathcal{R}}(\boldsymbol{U}, p), \quad \mathrm{where} \ {\mathcal{R}}(\boldsymbol{U}, p) = {\mathcal{D}}(\boldsymbol{U}) + \langle\!\langle p , \operatorname{DIV} \boldsymbol{U} \rangle\!\rangle. \end{equation}On a closed surface, the stationary conditions yield the Stokes equations for an evolving surface:
\begin{equation} 2 \mu \operatorname{DIV}^{\boldsymbol{\nabla}} {\mathcal{K}} \boldsymbol{U} - \operatorname{GRAD} p + \boldsymbol{B} - \delta_{\varphi} V_\varphi = 0,\quad \operatorname{DIV} \boldsymbol{U} = 0. \end{equation}
The symmetric, negative-definite, Laplacian-like operator 
$2\operatorname {DIV}^{\boldsymbol {\nabla }} {\mathcal {K}}: \varGamma (TW) \rightarrow \varGamma (T^*W)$ quantifies the viscous force, and is thus termed viscosity Laplacian. Although not the main focus of this paper, a common free energy in this context is the Helfrich Hamiltonian 
$V_{\varphi } = \int _W \kappa H^2\, \mathrm {d} A$, which results in the bending force 
$\boldsymbol {B}_\kappa = -\delta _\varphi V_\varphi = \kappa \boldsymbol {N} ( \Delta H + 2 H ( H^2 - K) )$, where H and K are the mean and Gaussian curvatures, respectively, 
$\kappa$ is the bending modulus, and 
$\boldsymbol {N}$ is the surface normal (Helfrich Reference Helfrich1973). This bending-driven relaxation dynamics is also used in the numerical examples of § 4. The commonly used form of (2.10), expressed through normal–tangent decomposition, is derived in the Appendix.
2.3.1. Remarks on the evolving Stokes equations and their solvability
Note that the presence of a ‘Killing vector field’ 
$\boldsymbol {X}$ (
${\mathcal {K}} \boldsymbol {X} = {{\boldsymbol{\mathsf{0}}}}$) can make the solution to (2.10) non-unique. In that case we select the least-norm solution, effectively projecting out rigid motions. This treatment aligns with scenarios where the surface is immersed with friction in a bulk fluid. A simple way to model friction is by augmenting the Rayleighian to 
${\mathcal {R}}^\alpha = {\mathcal {R}} + \alpha \| \boldsymbol {U} - \boldsymbol {U}_0 \|^2 / 2$, where 
$\boldsymbol {U}_0$ is the bulk velocity and 
$\alpha$ is the friction coefficient. This introduces a friction force 
$\boldsymbol {B}_\alpha = \alpha (\boldsymbol {U} - \boldsymbol {U}_0)$ in the momentum equation. Such flow driven by the bulk fluid is used as a numerical example in § 4.3. When 
$\boldsymbol {U}_0 = \boldsymbol {0}$, the solution approaches the least-norm solution as 
$\alpha \rightarrow 0$.
3. Methods
In this section, we present a simple numerical method to solve (2.10) on a triangle mesh with general geometry and topology. We begin by discretizing the strain rate tensor using finite elements, following the definition provided in (2.2). Based on the discretized strain rate, we construct the system's Rayleighian. A nonlinearly stable, structure-preserving variational integrator is derived based on Onsager's variational principle. To keep the numerical scheme minimal, we demonstrate the method using first-order spatial and temporal discretizations.
We discretize the system on a closed triangular mesh 
$M$ of arbitrary genus, consisting of vertices, edges and faces 
$\{\mathfrak {V}, \mathfrak {E}, \mathfrak {F}\}$. The vertex positions are given by the realization 
$\boldsymbol {\varphi }: M \hookrightarrow \mathbb {R}^3$. Here 
$M$ has a tangent space 
$T_\alpha M$ and a surface normal 
$\boldsymbol {N}_\alpha$ at each face 
$\alpha \in \mathfrak F$ (cf. figure 1).
Figure 1. Discretization of the strain rate tensor. The velocity 
$\boldsymbol {U}_i$ is defined at each vertex 
$i$, with surface normals 
$\boldsymbol {N}_\alpha$ at each face 
$\alpha$. The finite-element hat function 
$\varPhi _j$, shown by the colourmap, has local support around vertex 
$j$. Differentiating 
$\boldsymbol {U}_i$ via 
$\varPhi _j$ produces 
$(\boldsymbol {\nabla } \boldsymbol {U})_\beta$, which is symmetrized and projected using 
${\mathcal {P}}_\beta = {{\boldsymbol{\mathsf{g}}}} - \boldsymbol {N}_\beta \otimes \boldsymbol {N}_\beta$ to yield the strain rate tensor 
${{\boldsymbol{\mathsf{E}}}}_\beta$ on each face 
$\beta$.
In the discrete setting, we express tensor-valued functions and operators using index notation, with upper indices denoting their components in Cartesian coordinates. Vertex and face elements are denoted with lower indices, where vertex indices are Roman letters and face indices are Greek letters. For example, we denote vertex positions as 
$\varphi _j^i$, face normals as 
$N_\alpha ^i$, face pressure as 
$p_\alpha$, for 
$j \in \mathfrak V$, 
$\alpha \in \mathfrak F$ and 
$i \in \{1, 2, 3\}$.
3.1. Strain rate tensor and differential operators
To express the strain rate tensor 
${{\boldsymbol{\mathsf{E}}}}$ on 
$M$ (cf. (2.2)) in Cartesian coordinates, we first realize it in 
$\mathbb {R}^3$. We denote this realization as 
$2 \tilde{{{\boldsymbol{\mathsf{E}}}}} = {\mathcal {P}} (({\boldsymbol {\nabla }} \boldsymbol {U})^\top {{\boldsymbol{\mathsf{g}}}} + {{\boldsymbol{\mathsf{g}}}} {\boldsymbol {\nabla }} \boldsymbol {U} ) {\mathcal {P}}$, which satisfies 
$\tilde{{{\boldsymbol{\mathsf{E}}}}} [\![ \boldsymbol {V}, \boldsymbol {W} ]\!] = {{\boldsymbol{\mathsf{E}}}} [\![ {\mathcal {P}} \boldsymbol {V}, {\mathcal {P}} \boldsymbol {W} ]\!]$ for 
$\boldsymbol {V}, \boldsymbol {W} \in \varGamma (\mathbb {R}^3)$. Here, the linear projector 
${\mathcal {P}} = {{\boldsymbol{\mathsf{g}}}} - \boldsymbol {N} \otimes \boldsymbol {N}: \varGamma (\mathbb {R}^3) \rightarrow \varGamma (TM)$ projects 
$\mathbb {R}^3$ vectors onto the surface. In the discrete setting (cf. figure 1), the strain rate tensor 
$\tilde E^{kj}_\alpha$ and Killing operator 
$\mathcal {K}_{\alpha l}^{kj i}$ are expressed in terms of the velocity 
$U_{j}^{i}$, projector 
${\mathcal {P}}_{\alpha }^{ij} = \delta _{\alpha }^{ij} - N_{\alpha }^i N_{\alpha }^j$, and 
$\mathbb {R}^3$ component-wise surface derivative 
$({\boldsymbol {\nabla }}_{{\mathcal {P}}})_{\alpha j}^i$ as
\begin{equation} 2 \tilde E^{kj}_\alpha = \sum_{l, i} 2 {\mathcal{K}}_{\alpha l}^{kj i} U_{l}^{i} = \sum_{l, i} ( {\mathcal{P}}_{\alpha}^{ik} ({\boldsymbol{\nabla}}_{{\mathcal{P}}})^j_{\alpha l} + {\mathcal{P}}_{\alpha}^{ij} ({\boldsymbol{\nabla}}_{{\mathcal{P}}})^k_{\alpha l} )U_{l}^i. \end{equation}
Component-wise, 
$({\boldsymbol {\nabla }}_{{\mathcal {P}}})_{\alpha j}^i = {{\boldsymbol {\nabla }}^i \varPhi _{\alpha j} }$ follows scalar finite-element discretization, where 
$\varPhi _{\alpha j}$ is the hat function at vertex 
$j$ restricted to face 
$\alpha$. The divergence (2.4) is given by the trace 
$\operatorname {DIV}_{\alpha l}^i = \sum _j {\mathcal {K}}^{jj i}_{\alpha l}$. The gradient (2.5) and covariant divergence (2.6) can be obtained, up to a sign flip, by permuting the lower indices of the divergence and Killing operators, respectively.
3.2. Variational integrator by Onsager's variational principle
Given the vertex position 
$\boldsymbol {\varphi }_{(0)}$ and pressure 
$p_{(0)}$ at the current time 
$t = 0$, with a time step 
$\epsilon > 0$, the velocity is discretized as 
$\boldsymbol{U} = \mathring{\boldsymbol{\varphi }}/\epsilon = (\boldsymbol{\varphi } - \boldsymbol{\varphi }_{(0)})/\epsilon$, and the power input as 
$\dot V_{\boldsymbol {\varphi }} = (V_{\boldsymbol {\varphi }} - V_{\boldsymbol {\varphi }_{(0)}}) / \epsilon$. State variables 
$\boldsymbol {\varphi }_{(\epsilon )}$ and 
$p_{(\epsilon )}$ at 
$t = \epsilon$ follow a time-incremental Onsager's variational principle on a discrete Rayleighian (cf. § 2.3), 
$\boldsymbol {\varphi }_{(\epsilon )}, p_{(\epsilon )} = \arg \min _{\boldsymbol {\varphi }} \max _p {\mathcal {R}}$, where
\begin{align}
{\mathcal{R}}(\boldsymbol{\varphi}, p;
\boldsymbol{\varphi}_{(0)}, p_{(0)}, \epsilon) &:=
\epsilon^2 [ \mu \| \tilde{{{\boldsymbol{\mathsf{E}}}}}
\|^2 + \langle\!\langle \operatorname{DIV} \boldsymbol{U},
p \rangle\!\rangle + \dot V_{\boldsymbol{\varphi}} -
\langle\!\langle \boldsymbol{U}\,|\, \boldsymbol{B}
\rangle\!\rangle ]\nonumber\\ &\equiv \mu
{\mathcal{K}}_{\lambda m}^{kji}(A_{\lambda \gamma})^{{-}1}
{\mathcal{K}}_{ \gamma n }^{kj l} \mathring{\varphi}|^i_m
\mathring{\varphi}|_n^{l} + \epsilon (\mathring{\varphi}|_{i}^l \operatorname{DIV}_{\alpha i}^l p_{\alpha} +
V_{\boldsymbol{\varphi}} - \mathring \varphi|_{i}^l B_{i}^l
)\nonumber\\ &\equiv \mu \mathring{\boldsymbol{\varphi}}^\top {{\boldsymbol{\mathsf{L}}}}
~\mathring{\boldsymbol{\varphi}} + \epsilon [
(\operatorname{DIV} \mathring{\boldsymbol{\varphi}})^\top p
+ V_{\boldsymbol{\varphi}} - \mathring{\boldsymbol{\varphi}}^\top \boldsymbol{B} ].
\end{align}
Here, we assume Einstein summation and 
${{\boldsymbol{\mathsf{L}}}} = {\mathcal {K}}^\top {{\boldsymbol{\mathsf{A}}}}^{-1} {\mathcal {K}}$ is half of the discrete viscosity Laplacian (cf. § 2.3), where 
$A_{\lambda \gamma }$ is the area of the face 
$\gamma$ for 
$\lambda =\gamma$ and zero otherwise.
Under any choice of 
$V_{\boldsymbol {\varphi }}$ (assuming 
$\boldsymbol {B} = \boldsymbol 0$), the variational integrator is unconditionally stable by design, as it preserves the system's dissipative structure. The dynamics follow a gradient flow of 
$V_{\boldsymbol {\varphi }}$ in the metric space defined by the Stokes operator, ensuring a monotonic decrease in the Rayleighian: 
${\mathcal {R}}_{(\epsilon )} \leq {\mathcal {R}}_{(0)}$. Therefore, the allowable time step size depends only on the solvability of the optimization in (3.2), which can be efficiently handled using standard numerical optimization methods. Here we use a simple gradient flow to solve this saddle-point optimization: 
${{\boldsymbol{\mathsf{H}}}} \dot{\boldsymbol {\varphi }} = - \delta _{\boldsymbol {\varphi }} {\mathcal {R}} = - 2 \mu { {{\boldsymbol{\mathsf{L}}}} \mathring{\boldsymbol {\varphi }} } - \epsilon ({\operatorname {DIV}^\top p } + {\delta _{\boldsymbol {\varphi }} V_{\boldsymbol {\varphi }}} - {\boldsymbol {B}})$ and 
$h \dot p = \delta _{p} \mathcal {R} = \operatorname {DIV} \mathring{\boldsymbol {\varphi }}$, where 
${{\boldsymbol{\mathsf{H}}}}$ and 
$h$ provide a general metric for the gradient flow. Remeshing is performed using SideFx Houdini.) In § 4, we adopt the Helfrich Hamiltonian as 
$V$ and follow the discretization used in Zhu, Lee & Rangamani (Reference Zhu, Lee and Rangamani2022).
As mentioned in § 2.3.1, (2.10) and its discrete version (3.2) have non-unique solutions up to rigid body modes in 
$\mathbb {R}^3$. To ensure that we obtain continuous dynamics during time evolution, we consistently project out rigid body modes.
A MATLAB implementation of our discrete theory is available at https://t.ly/vUxfh. Tensorial calculations in (3.1) and (3.2) utilize 
$\texttt{sptensor}$ (Bader & Kolda Reference Bader and Kolda2008), with remeshing performed using SideFx Houdini as needed.
4. Results
In this section, we validate the discrete model, analyse its convergence properties, and demonstrate its applicability to manifolds with intricate geometries and topologies. Numerical validation against analytical solutions is presented in § 4.1, with further examples provided in §§ 4.2 and 4.3 to showcase its effectiveness.
4.1. Validations
We validate the differential operators 
$\mathcal{K}$ and 
$\operatorname{GRAD}$ on a spheroid with semi-axes a and c, parametrized by latitude 
$\beta \in [- {\rm \pi}/2, {\rm \pi}/2]$ and longitude 
$\theta \in (-{\rm \pi}, {\rm \pi}]$, using the realization 
$\boldsymbol{\varphi}(\beta, \theta) = [a \cos \beta \cos \theta, a \cos \beta \sin \theta, c \sin \beta]$.
In figure 2(a), we test the operator 
$\mathcal {K}$ by computing the lowest 10 eigenvalues of 
${{\boldsymbol{\mathsf{L}}}} = \mathcal {K}^\top {{\boldsymbol{\mathsf{A}}}}^{-1}{\mathcal {K}}$ with 
${{\boldsymbol{\mathsf{L}}}} \boldsymbol {U}_i = \lambda _i {{\boldsymbol{\mathsf{A}}}} \boldsymbol {U}_i$ (cf. (3.2)). As predicted by the continuous theory, there are six vanishing eigenvalues, forming a six dimensional kernel of 
${\mathcal {K}}$ corresponding to the space of rigid body motions in 
$\mathbb {R}^3$.
Figure 2. Validation of our model and method. (a) Minimum eigenvalues 
$\lambda _{1\leq i \leq 10}$ of 
${{\boldsymbol{\mathsf{L}}}} \boldsymbol {U}_i = \lambda _i {{\boldsymbol{\mathsf{A}}}} \boldsymbol {U}_i$. (b) Comparison of numerical 
$\operatorname {GRAD} 1$ and analytical 
$H(\beta )$ evaluations of the mean curvature on a spheroid. (c) Relative error 
$\varepsilon _{\boldsymbol {U}}$ for the momentum equation and 
$\varepsilon _{p}$ for the continuity equation as functions of mean edge length 
$l$, with the inset showing the tangential forcing 
$\boldsymbol {b}$ constructed from a spherical harmonic. (d) Relative error in total area A compared to the initial area 
$A_0$ during Helfrich–Stokes relaxation, with insets showing the initial frame (prolate spheroid) and the final frame (sphere) of the relaxation.
To validate the operator 
$\operatorname {GRAD}$ (A4), we check that the numerical evaluation of 
${{\boldsymbol{\mathsf{A}}}}^{-1}|\operatorname {DIV}^\top 1|$ aligns with the analytical mean curvature 
$H(\beta )$ of an oblate spheroid with 
$a = 1$ and 
$c = 0.5$, as shown in figure 2(b). In fact, this evaluation of mean curvature is identical with the established cotan Laplacian discretization (Meyer et al. Reference Meyer, Desbrun, Schröder and Barr2003).
We evaluate the convergence of our discretization of (2.10) by comparing it to an analytical solution on a unit sphere. The forcing term is prescribed analytically as 
$\boldsymbol {B} = \boldsymbol {b} + b_n \boldsymbol {N} = \boldsymbol {N} \times \operatorname {grad} \phi + b_n \boldsymbol {N}$, where 
$\phi$ is a spherical harmonic satisfying 
$\Delta \phi := \operatorname {div} \operatorname {grad} \phi = \lambda \phi$, 
$\boldsymbol {N}$ is the surface normal, and 
$b_n$ is a constant. On the unit sphere, 
$\boldsymbol {b}$ is an eigenfunction of the viscosity Laplacian with eigenvalue 
$\lambda + 2$. The analytical solution 
$\overline {\boldsymbol {U}} = -\boldsymbol {b} / (\lambda + 2)$ and 
$\bar {p} = -b_n / 2$ satisfies (2.10) (
$\mu = 1$, 
$\delta _{\varphi } V_{\varphi } = 0$) exactly (cf. (A7)). Numerically, the relative residual of the momentum equation, 
$\varepsilon _{\boldsymbol {U}} = \|2 {{\boldsymbol{\mathsf{A}}}}^{-1} {{\boldsymbol{\mathsf{L}}}} \overline {\boldsymbol {U}} - {{\boldsymbol{\mathsf{A}}}}^{-1} \operatorname {DIV}^\top \bar {p} - \boldsymbol {B} \| / \|\boldsymbol {B} \|$, and the residual of the continuity equation, 
$\varepsilon _p = \|{{\boldsymbol{\mathsf{A}}}}^{-1} \operatorname {DIV} \overline{\boldsymbol {U}} \|$, shown in figure 2(c), decrease approximately linearly as the mesh is refined.
Lastly, figure 2(d) illustrates the total area of the fluid film during the Stokes relaxation (
$\mu = 1$) from a prolate spheroid (
$a = 1$, 
$c = 2$) to a sphere under Helfrich energy (
$\kappa = 0.05$), without mesh reparametrization. Local incompressibility conserves the global area within a relative error of 
$10^{-5}$ throughout the evolution.
4.2. Relaxation of a genus-6 torus
We model the evolution of a high-genus, non-analytical lipid membrane (
$\mu = 1$, 
$\kappa = 0.05$) based on the Helfrich–Stokes relaxation in figure 3(a). In theory, the elastic Helfrich Hamiltonian should be monotonically dissipated through viscosity, 
$V + E_{\mu } = \int \kappa H^2 \,\mathrm {d}A + \int _0^t \langle \!\langle {{\boldsymbol{\mathsf{T}}}} \,|\, {{\boldsymbol{\mathsf{E}}}} \rangle \!\rangle / 2|_\tau \,\mathrm {d}\tau = \mathrm {const}$. Although our numerical results in figure 3(b) do not exactly match this theoretical expectation, they show very good agreement.
Figure 3. (a) Helfrich–Stokes relaxation of a genus-6 torus. Snapshots are taken at 
$t = 0$, 
$2$, 
$8$, 
$14$, 
$20$, and 
$40$, from left to right, top to bottom. The animated simulation is available at https://youtu.be/Llh0_N0hCPw and in movie 1 of the supplementary movies available at https://doi.org/10.1017/jfm.2024.1208. (b) Elastic energy 
$V = \int \kappa H^2 \,\mathrm {d}A$ and cumulative dissipation 
$E_{\mu } = \int _0^t \langle \!\langle {{\boldsymbol{\mathsf{T}}}} \,|\, {{\boldsymbol{\mathsf{E}}}} \rangle \!\rangle / 2|_\tau \,\mathrm {d}\tau$, which theoretically sum to a constant value.
4.3. Deformation of a lipid membrane under bulk flow
We model the deformation of a spherical lipid vesicle (
$\mu = 1$, 
$\kappa = 0.01$) under the friction (
$\alpha = 0.5$, cf. § 2.3.1) from a constant bulk shear-extension flow defined by 
$\boldsymbol {U}_0(\theta, r, z) = \boldsymbol {U}_0^{shear} + \boldsymbol {U}_0^{ext} = (5 r z) \hat{\boldsymbol{e}}_\theta + (- 0.5 r \hat{\boldsymbol{e}}_r + z \hat{\boldsymbol{e}}_z)$ in cylindrical coordinates. As illustrated in figure 4(a), this bulk flow stretches and twists the initially spherical vesicle along the 
$z$-direction. A snapthrough-like instability occurs around 
$t = 3.7$ between frame (iii) and frame (iv) in figure 4(a). While the bulk friction governs most of the dynamics, the Helfrich energy dominates during this buckling event, as indicated by the dip of the energy dissipation rate in figure 4(b) around 
$t = 3.7$.
Figure 4. (a) Deformation of a lipid membrane under bulk flow, with the tangent velocity 
$\boldsymbol {u}$ shown as streamlines and the normal velocity 
$U_{\boldsymbol {N}}$ displayed in a colourmap ranging from 
$-1$ to 
$1$ (cf. the Appendix for the tangent–normal decomposition). Snapshots (i)–(iv) are taken at 
$t = 1.5$, 
$2.5$, 
$3.4$ and 
$3.8$, respectively. The animated simulation is available at https://youtu.be/M0WsLihzRJk and in movie 2 of the supplementary movies. (b) Energy dissipation rate 
$\dot E$ of the system due to viscosity 
$\dot {E}_{\mu } = \langle \!\langle {{\boldsymbol{\mathsf{T}}}}\,|\, {{\boldsymbol{\mathsf{E}}}} \rangle \!\rangle / 2$, elasticity 
$\dot {E}_{\kappa } = \langle \!\langle \boldsymbol {B}_\kappa \,|\, \boldsymbol {U} \rangle \!\rangle$ and bulk friction 
$\dot {E}_{\alpha } = \langle \!\langle \boldsymbol {B}_\alpha \,|\, \boldsymbol {U} \rangle \!\rangle / 2$.
5. Conclusion
Using Onsager's variational principle, we derived and reformulated the evolving Stokes equations on a manifold. We replaced the classic system of equations expressed in tangent–normal splitting with a simple, coordinate-free differential-geometric formulation. This approach leads directly to a straightforward discrete model and a numerical scheme to solve this long-standing problem in its full geometric generality.
Although the theoretical framework is not limited, the current minimal implementation is restricted to closed manifolds and employs first-order discretization. Future work could incorporate boundary conditions and higher-order discretizations. For applications, we expect the minimal system presented here to serve as a foundation for integrating more complex models for specific biophysical problems. These could include additional global volumetric/areal constraints, heterogeneity in material properties, in-plane anisotropy and surface activity, as exemplified in Zhu et al. (Reference Zhu, Lee and Rangamani2022) and Zhu, Saintillan & Chern (Reference Zhu, Saintillan and Chern2024).
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2024.1208.
Funding
The authors acknowledge funding from National Science Foundation grant no. DMS-2153520 (D.S.) and National Science Foundation CAREER Award 2239062 (A.C.). Additional support was provided by SideFX software.
Declaration of interests
The authors report no conflict of interest.
Appendix. Evolving Stokes equations through normal–tangent decomposition
A.1. Differential operators for an evolving surface
The second fundamental form (i.e. curvature tensor), 
$\boldsymbol{\mathsf{II}} \in \varGamma (T^*M \otimes _{symm}T^*M)$, is given by the derivative of the surface normal 
$\boldsymbol {N}$:
\begin{equation} \boldsymbol{\mathsf{II}} [\![ \boldsymbol{v}, \boldsymbol{w} ]\!] =( {{\boldsymbol{\mathsf{F}}}} \boldsymbol{v} )\boldsymbol{\cdot}({\boldsymbol{\nabla}}_{{{\boldsymbol{\mathsf{F}}}} \boldsymbol{w}} \boldsymbol{N}) =:\langle \boldsymbol{v}, {{\boldsymbol{\mathsf{S}}}} \boldsymbol{w} \rangle, \end{equation}
where 
${{\boldsymbol{\mathsf{S}}}}: \varGamma (TM) \rightarrow \varGamma (TM)$ is the shape operator.
With the induced metric 
$\boldsymbol{\mathsf{I}}$, there is an intrinsic Levi-Civita covariant derivative 
${\boldsymbol {\nabla }}: \varGamma (TM) \rightarrow \varGamma (TM \otimes T^*M)$ on 
$M$. The covariant derivative in 
$W$ used in (A1) and (2.1) can be expressed in terms of the intrinsic covariant derivative and curvature tensor as
\begin{equation} {\boldsymbol{\nabla}}_{{{\boldsymbol{\mathsf{F}}}} \boldsymbol{w}} {{\boldsymbol{\mathsf{F}}}} {\boldsymbol{v}} = {\boldsymbol{\nabla}}_{\boldsymbol{w}} \boldsymbol{v} - \boldsymbol{\mathsf{II}} [\![ \boldsymbol{w}, \boldsymbol{v} ]\!] \boldsymbol{N}. \end{equation} We can decompose the velocity into tangent and normal parts, 
$\boldsymbol {U} = {{\boldsymbol{\mathsf{F}}}} {\boldsymbol {u}} + U_N \boldsymbol {N}$, where 
$\boldsymbol {u} \in \varGamma (TM)$, 
$U_N \in C^\infty (W)$. Since 
${{\boldsymbol{\mathsf{F}}}}{\boldsymbol {u}}$ effectively embeds the tangent vector 
$\boldsymbol {u}$ in 
$\mathbb {R}^3$, we will abbreviate 
${{\boldsymbol{\mathsf{F}}}} {\boldsymbol {u}}$ as 
$\boldsymbol {u}$ when unambiguous. Given 
$\boldsymbol {v} \in T_pM$, (2.1) can be expressed as 
${\boldsymbol {\nabla }}_{{{\boldsymbol{\mathsf{F}}}} \boldsymbol {v}} \boldsymbol {U} = { {\boldsymbol {\nabla }}_{\boldsymbol {v}} \boldsymbol {u} + U_N {{\boldsymbol{\mathsf{S}}}} \boldsymbol {v}} + \boldsymbol {N} ( \langle \operatorname {grad} U_N, \boldsymbol {v} \rangle - \boldsymbol{\mathsf{II}} [\![ \boldsymbol {u}, \boldsymbol {v} ]\!] )$. Combining with (2.2) and (2.3), we have 
$2 ({\mathcal {K}} \boldsymbol {U}) [\![ \boldsymbol {v}, \boldsymbol {w} ]\!] = \langle {\boldsymbol {\nabla }}_{\boldsymbol {v}} \boldsymbol {u}, \boldsymbol {w} \rangle + \langle {\boldsymbol {\nabla }}_{\boldsymbol {w}} \boldsymbol {u} , \boldsymbol {v} \rangle + 2 U_N \boldsymbol{\mathsf{II}} [\![ \boldsymbol {v}, \boldsymbol {w} ]\!]$, or
\begin{equation} 2 {\mathcal{K}} \boldsymbol{U}= {\boldsymbol{\nabla}} \boldsymbol{u} + ({\boldsymbol{\nabla}} \boldsymbol{u})^\top + 2 U_N \boldsymbol{\mathsf{II}} =: 2 \mathscr{K} \boldsymbol{u} + 2 U_N \boldsymbol{II}, \end{equation}
where 
$\mathscr {K}: \varGamma (TM) \rightarrow \varGamma (T^*M \otimes _{\operatorname {symm}} T^*M)$ denotes the intrinsic Killing operator. Combining (2.4), (2.5) and (A3), the divergence and gradient can be expressed as
\begin{equation} \operatorname{DIV} \boldsymbol{U} = \operatorname{div}\boldsymbol{u} + 2H U_N , \quad \operatorname{GRAD}(p) = \operatorname{grad} p - 2 p H \boldsymbol{N}, \end{equation}
where 
$\operatorname {div} = \operatorname {tr} (\circ \mathscr {K})$ is the intrinsic divergence, 
$\operatorname {grad} = -\operatorname {div}^*$ is the intrinsic gradient, and 
$H = \operatorname {tr} \boldsymbol{\mathsf{II}} /2$ is the mean curvature of the surface. From (2.6) and (A3), we get:
\begin{equation} \operatorname{DIV}^{\boldsymbol{\nabla}} {{\boldsymbol{\mathsf{T}}}} = \operatorname{div}^{\boldsymbol{\nabla}} {{\boldsymbol{\mathsf{T}}}} - \langle \boldsymbol{\mathsf{II}} \,| \,{{\boldsymbol{\mathsf{T}}}} \rangle \boldsymbol{N}, \end{equation}
where 
$\operatorname {div}^{\boldsymbol {\nabla }} = -\mathscr K^*: \varGamma (TM \otimes TM) \rightarrow \varGamma (T^*M)$ is the intrinsic covariant divergence.
A.2. Evolving Stokes equations
The viscosity Laplacian for an evolving surface can be explicitly decomposed as
\begin{align} 2 \operatorname{DIV}^{\boldsymbol{\nabla}} {\mathcal{K}} \boldsymbol{U} &= 2 [ \operatorname{div}^{\boldsymbol{\nabla}} \mathscr{K} \boldsymbol{u} + \operatorname{div}^{\boldsymbol{\nabla}} (U_N\boldsymbol{\mathsf{II}}) - ( \langle \boldsymbol{\mathsf{II}} , {\mathscr{K} \boldsymbol{u}} \rangle + U_N |\boldsymbol{\mathsf{II}}|^2 ) \boldsymbol{N} ]\nonumber\\ &= ( \boldsymbol{\Delta} + K + \operatorname{grad} \circ \operatorname{div})\boldsymbol{u} + 2 {{\boldsymbol{\mathsf{S}}}} (\operatorname{grad} U_N) + 4 U_N \operatorname{grad} H\nonumber\\ &\quad -2 [ \langle {\boldsymbol{\nabla}} \boldsymbol{u}, \boldsymbol{\mathsf{II}} \rangle + U_N (4 H^2 - 2K) ] \boldsymbol{N}, \end{align}
where 
$\boldsymbol {\Delta } = \operatorname {div}^{\boldsymbol {\nabla }} {\boldsymbol {\nabla }}: \varGamma (TM) \rightarrow \varGamma (T^*M)$ is the intrinsic Bochner Laplacian, and 
$K = \det \boldsymbol{\mathsf{II}}$ is the Gaussian curvature.
To summarize, by substituting (A3)–(A6) into (2.10), we get the system of evolving Stokes equations subject to external force 
$\boldsymbol {B} = \boldsymbol {b} + b_n \boldsymbol {N}$ and bending force 
$\boldsymbol {B}_\kappa = \kappa \boldsymbol {N} ( \Delta H + 2 H ( H^2 - K) )$:
\begin{align} \left.\begin{array}{c} \mu
[ (\boldsymbol{\Delta} + K)\boldsymbol{u} + 2
{{\boldsymbol{\mathsf{Q}}}} (\operatorname{grad} U_N) + 2 U_N
\operatorname{grad} H ] - \operatorname{grad} p +
\boldsymbol{b} = 0, \\ -2 \mu [ \langle
{\boldsymbol{\nabla}} \boldsymbol{u} , \boldsymbol{\mathsf{II}}
\rangle + U_N (4H^2 - 2K) ] + 2 p H + b_n + \kappa ( \Delta
H + 2 H ( H^2 - K) ) = 0,\\ \operatorname{div}
\boldsymbol{u} + 2H U_N = 0. \end{array}\right\}
\end{align}
Note that the Hopf differential 
${{\boldsymbol{\mathsf{Q}}}} = {{\boldsymbol{\mathsf{S}}}} - H \boldsymbol{\mathsf{I}}$ is the traceless part of the shape operator 
${{\boldsymbol{\mathsf{S}}}}$ that represents the deviatoric curvature. Formulation (A7) agrees with earlier results by Scriven (Reference Scriven1960) and Arroyo & DeSimone (Reference Arroyo and DeSimone2009).
Open access
