Approximate Mean Curvature Flow

• Approximate Mean Curvature Flow is a class of methods that regularize and discretize the classical curvature flow to handle singular, degenerate, or nonsmooth data.
• Key approaches include level-set regularization, diffusion-thresholding, diffuse-interface limits, variational finite element schemes, and varifold-based approximations.
• These techniques offer practical insights into controlling discretization and regularization errors, ensuring numerical convergence and robustness in complex geometric evolutions.

Approximate mean curvature flow denotes a family of regularized, discretized, and weakly formulated evolutions that are designed to represent motion by mean curvature when the classical flow is singular, degenerate, or defined on nonsmooth data. In the literature, approximation appears through elliptic level-set regularization, heat-kernel thresholding, diffuse-interface limits, stochastic control representations, geometric finite element schemes, point-cloud and varifold dynamics, and smooth surgery constructions. These approaches target different notions of solution—viscosity solutions, semigroup limits, energy-dissipative discrete flows, Brakke motions, or spacetime varifold limits—and they also extend to power-law curvature speeds, obstacle constraints, sub-Riemannian settings, and boundary value problems (Kröner et al., 2015, Kröner, 2013, Baspinar et al., 2016, Pisante et al., 2013, Buet et al., 1 Oct 2025).

1. Formulations and target notions of convergence

Approximation frameworks for mean curvature flow are organized by the representation of the interface and by the limiting notion of motion. In level-set methods, the evolving hypersurface is encoded as a level set of a scalar function and the singular geometric operator is regularized before discretization. In diffusion-thresholding methods, short-time heat evolution is alternated with a geometric reconstruction step. In diffuse-interface methods, a reaction-diffusion PDE generates energy measures whose sharp-interface limit is a Brakke flow. In geometric finite element methods, the curve or surface itself is discretized and evolved by a structure-preserving ODE or DAE system. In varifold-based approaches, approximate mean curvature is defined directly from the first variation and mass of a measure-theoretic surface, allowing arbitrary dimension, codimension, and even discrete point clouds (Kröner et al., 2015, Baspinar et al., 2016, Pisante et al., 2013, Liu et al., 2024, Buet et al., 1 Oct 2025).

Paradigm	Representative mechanism	Typical limit notion
Level-set regularization	Elliptic regularization of singular curvature PDE	Viscosity solution / FE approximation
Diffusion-thresholding	Heat equation plus thresholding or directional reconstruction	Nonlinear semigroup / motion by curvature
Diffuse-interface	Allen-Cahn or related phase-field evolution	Brakke motion
Geometric discretization	Variational or FE evolution of nodes or meshes	Discrete energy-dissipative flow
Varifold approximation	Regularized first variation and push-forward dynamics	Brakke-type equality / spacetime Brakke flow

A central point is that “approximate mean curvature flow” is not a single model. It is a methodological class whose members preserve different structural features: consistency with viscosity theory, monotonicity of mass or energy, weak compactness, geometric comparison, or implementability on irregular data. This suggests a taxonomy by solution concept rather than by a single PDE.

2. Level-set regularization and finite element approximation

A widely used route starts from a level-set formulation. For power mean curvature flow with normal speed $H^k$, $k \ge 1$, the time-dependent level-set equation is

$$\frac{\partial u}{\partial t} = |Du| \left( \operatorname{div}\frac{Du}{|Du|} \right)^k.$$

For $k>1$, the time-dependent formulation presents analytical difficulties, and the stationary level-set formulation attributed to Schulze is used: $$\operatorname{div}\left(\frac{Du}{|Du|}\right) = -\frac{1}{|Du|^{1/k}} \quad \text{in } \Omega, \qquad u=0 \quad \text{on } \partial\Omega.$$ Its elliptic regularization is

$$\operatorname{div}\left(\frac{Du^\varepsilon}{\sqrt{\varepsilon^2+|Du^\varepsilon|^2}}\right) = -(\varepsilon^2+|Du^\varepsilon|^2)^{-1/(2k)} \quad \text{in } \Omega, \qquad u^\varepsilon=0 \quad \text{on } \partial\Omega,$$

which removes the singularity at $Du=0$ (Kröner et al., 2015).

Finite element discretization is then performed on a triangulated or tetrahedral mesh. One formulation uses continuous linear elements

$$V_h=\{v\in C^0(\bar\Omega_h): v|_{\partial\Omega_h}=0,\ v|_T \text{ linear}\},$$

and seeks $u_h^\varepsilon\in V_h$ satisfying the nonlinear weak problem

$$\int_{\Omega_h} \frac{Du_h^\varepsilon\cdot Dv_h}{\sqrt{\varepsilon^2+|Du_h^\varepsilon|^2}}\,dx = \int_{\Omega_h} (\varepsilon^2+|Du_h^\varepsilon|^2)^{-1/(2k)} v_h\,dx \qquad \forall v_h\in V_h.$$

A related earlier analysis for $k \ge 1$0 used continuous finite elements that are polynomials of degree $k \ge 1$1 on each tetrahedron, proved existence and uniqueness of $k \ge 1$2 sufficiently close to $k \ge 1$3, and indicated that for $k \ge 1$4 higher-order finite elements are required (Kröner, 2013).

The approximation error is split into regularization and discretization parts: $k \ge 1$5 For fixed $k \ge 1$6, one obtains

$k \ge 1$7

while the regularization error satisfies, for $k \ge 1$8,

$k \ge 1$9

and more generally

$$\frac{\partial u}{\partial t} = |Du| \left( \operatorname{div}\frac{Du}{|Du|} \right)^k.$$0

The total error can therefore be controlled by coupling $$\frac{\partial u}{\partial t} = |Du| \left( \operatorname{div}\frac{Du}{|Du|} \right)^k.$$1 and $$\frac{\partial u}{\partial t} = |Du| \left( \operatorname{div}\frac{Du}{|Du|} \right)^k.$$2 polynomially (Kröner et al., 2015).

The numerical studies in the curve case confirm these rates and supply geometric information. For circular data in the $$\frac{\partial u}{\partial t} = |Du| \left( \operatorname{div}\frac{Du}{|Du|} \right)^k.$$3 case, the exact solution is

$$\frac{\partial u}{\partial t} = |Du| \left( \operatorname{div}\frac{Du}{|Du|} \right)^k.$$4

The experiments reported convergence in $$\frac{\partial u}{\partial t} = |Du| \left( \operatorname{div}\frac{Du}{|Du|} \right)^k.$$5, $$\frac{\partial u}{\partial t} = |Du| \left( \operatorname{div}\frac{Du}{|Du|} \right)^k.$$6, and $$\frac{\partial u}{\partial t} = |Du| \left( \operatorname{div}\frac{Du}{|Du|} \right)^k.$$7 norms, sometimes faster than worst-case estimates, and showed that the isoperimetric deficit $$\frac{\partial u}{\partial t} = |Du| \left( \operatorname{div}\frac{Du}{|Du|} \right)^k.$$8 decreases monotonically under the evolution. They also reported that as $$\frac{\partial u}{\partial t} = |Du| \left( \operatorname{div}\frac{Du}{|Du|} \right)^k.$$9 increases, the curve becomes round and shrinks more quickly to a point. Implementation used linear finite elements on unstructured meshes created with Gmsh and Newton-type solvers for the nonlinear systems (Kröner et al., 2015).

3. Diffusion, thresholding, stochastic control, and diffuse-interface limits

The Bence-Merriman-Osher algorithm approximates mean curvature flow by iterating short-time heat diffusion and thresholding. Starting from a set $k>1$0, one solves

$k>1$1

for a short time $k>1$2, then defines

$k>1$3

and repeats. Evans’ convergence argument reformulates the iteration through a nonlinear semigroup and level-set theory. A modified version replaces the initial indicator by the surface measure datum

$k>1$4

evolves it by the heat equation, and defines the new surface by

$k>1$5

For small $k>1$6, the induced normal velocity is

$k>1$7

and when $k>1$8 this yields the graph formulation of classical mean curvature flow up to $k>1$9. Global convergence is again proved by nonlinear semigroup methods (Baspinar et al., 2016).

A probabilistic approximation is available in sub-Riemannian geometry through a Riemannian approximation. The approximated mean curvature PDE is written as

$$\operatorname{div}\left(\frac{Du}{|Du|}\right) = -\frac{1}{|Du|^{1/k}} \quad \text{in } \Omega, \qquad u=0 \quad \text{on } \partial\Omega.$$0

or equivalently $$\operatorname{div}\left(\frac{Du}{|Du|}\right) = -\frac{1}{|Du|^{1/k}} \quad \text{in } \Omega, \qquad u=0 \quad \text{on } \partial\Omega.$$1. The associated controlled diffusion is

$$\operatorname{div}\left(\frac{Du}{|Du|}\right) = -\frac{1}{|Du|^{1/k}} \quad \text{in } \Omega, \qquad u=0 \quad \text{on } \partial\Omega.$$2

with value function

$$\operatorname{div}\left(\frac{Du}{|Du|}\right) = -\frac{1}{|Du|^{1/k}} \quad \text{in } \Omega, \qquad u=0 \quad \text{on } \partial\Omega.$$3

For fixed $$\operatorname{div}\left(\frac{Du}{|Du|}\right) = -\frac{1}{|Du|^{1/k}} \quad \text{in } \Omega, \qquad u=0 \quad \text{on } \partial\Omega.$$4, this value function is proved to be a bounded lower semicontinuous viscosity solution of the approximated mean curvature flow, and uniqueness follows under the comparison principle (Grande, 2021).

Diffuse-interface approximation proceeds instead through the Allen-Cahn equation on a Riemannian manifold: $$\operatorname{div}\left(\frac{Du}{|Du|}\right) = -\frac{1}{|Du|^{1/k}} \quad \text{in } \Omega, \qquad u=0 \quad \text{on } \partial\Omega.$$5 The associated energy density measures

$$\operatorname{div}\left(\frac{Du}{|Du|}\right) = -\frac{1}{|Du|^{1/k}} \quad \text{in } \Omega, \qquad u=0 \quad \text{on } \partial\Omega.$$6

converge, along subsequences, to rectifiable Radon measures evolving by mean curvature flow in the sense of Brakke. A local almost monotonicity formula, described as a weak counterpart of Huisken’s monotonicity formula, supplies density control and is a key tool in passing from the diffuse interface to the Brakke limit on space forms (Pisante et al., 2013).

Diffusion-thresholding also extends to constrained evolutions. For obstacle mean curvature flow, one augments the MBO step by a pointwise clamping rule: $$\operatorname{div}\left(\frac{Du}{|Du|}\right) = -\frac{1}{|Du|^{1/k}} \quad \text{in } \Omega, \qquad u=0 \quad \text{on } \partial\Omega.$$7 equivalently thresholding followed by the update $$\operatorname{div}\left(\frac{Du}{|Du|}\right) = -\frac{1}{|Du|^{1/k}} \quad \text{in } \Omega, \qquad u=0 \quad \text{on } \partial\Omega.$$8 on $$\operatorname{div}\left(\frac{Du}{|Du|}\right) = -\frac{1}{|Du|^{1/k}} \quad \text{in } \Omega, \qquad u=0 \quad \text{on } \partial\Omega.$$9, $$\operatorname{div}\left(\frac{Du^\varepsilon}{\sqrt{\varepsilon^2+|Du^\varepsilon|^2}}\right) = -(\varepsilon^2+|Du^\varepsilon|^2)^{-1/(2k)} \quad \text{in } \Omega, \qquad u^\varepsilon=0 \quad \text{on } \partial\Omega,$$0 on $$\operatorname{div}\left(\frac{Du^\varepsilon}{\sqrt{\varepsilon^2+|Du^\varepsilon|^2}}\right) = -(\varepsilon^2+|Du^\varepsilon|^2)^{-1/(2k)} \quad \text{in } \Omega, \qquad u^\varepsilon=0 \quad \text{on } \partial\Omega,$$1. The resulting scheme retains a geometric comparison principle and a minimizing movements interpretation, is unconditionally stable, and converges to the viscosity solution of obstacle mean curvature flow (Krämer et al., 18 Dec 2025).

4. Variational and geometric finite element discretizations

A distinct approximation philosophy keeps the interface explicit and discretizes the underlying gradient structure. Using the Onsager principle, mean curvature flow of a smooth closed curve in $$\operatorname{div}\left(\frac{Du^\varepsilon}{\sqrt{\varepsilon^2+|Du^\varepsilon|^2}}\right) = -(\varepsilon^2+|Du^\varepsilon|^2)^{-1/(2k)} \quad \text{in } \Omega, \qquad u^\varepsilon=0 \quad \text{on } \partial\Omega,$$2 is derived from the Rayleighian built from the length energy

$$\operatorname{div}\left(\frac{Du^\varepsilon}{\sqrt{\varepsilon^2+|Du^\varepsilon|^2}}\right) = -(\varepsilon^2+|Du^\varepsilon|^2)^{-1/(2k)} \quad \text{in } \Omega, \qquad u^\varepsilon=0 \quad \text{on } \partial\Omega,$$3

and the dissipation

$$\operatorname{div}\left(\frac{Du^\varepsilon}{\sqrt{\varepsilon^2+|Du^\varepsilon|^2}}\right) = -(\varepsilon^2+|Du^\varepsilon|^2)^{-1/(2k)} \quad \text{in } \Omega, \qquad u^\varepsilon=0 \quad \text{on } \partial\Omega,$$4

Minimization gives $$\operatorname{div}\left(\frac{Du^\varepsilon}{\sqrt{\varepsilon^2+|Du^\varepsilon|^2}}\right) = -(\varepsilon^2+|Du^\varepsilon|^2)^{-1/(2k)} \quad \text{in } \Omega, \qquad u^\varepsilon=0 \quad \text{on } \partial\Omega,$$5, with energy dissipation

$$\operatorname{div}\left(\frac{Du^\varepsilon}{\sqrt{\varepsilon^2+|Du^\varepsilon|^2}}\right) = -(\varepsilon^2+|Du^\varepsilon|^2)^{-1/(2k)} \quad \text{in } \Omega, \qquad u^\varepsilon=0 \quad \text{on } \partial\Omega,$$6

After polygonal discretization, the discrete flow becomes a system of ODEs

$$\operatorname{div}\left(\frac{Du^\varepsilon}{\sqrt{\varepsilon^2+|Du^\varepsilon|^2}}\right) = -(\varepsilon^2+|Du^\varepsilon|^2)^{-1/(2k)} \quad \text{in } \Omega, \qquad u^\varepsilon=0 \quad \text{on } \partial\Omega,$$7

and the semidiscrete energy law

$$\operatorname{div}\left(\frac{Du^\varepsilon}{\sqrt{\varepsilon^2+|Du^\varepsilon|^2}}\right) = -(\varepsilon^2+|Du^\varepsilon|^2)^{-1/(2k)} \quad \text{in } \Omega, \qquad u^\varepsilon=0 \quad \text{on } \partial\Omega,$$8

is preserved exactly. The method is advanced in time by the improved Euler scheme, extended to volume-preserving mean curvature flow and wetting problems, and numerical examples show optimal convergence rate for all three problems (Liu et al., 2024).

For open curves in two-dimensional conformally flat Riemannian manifolds, variational weak formulations lead to piecewise linear finite element schemes for curvature flow and elastic flow. The geodesic curvature is

$$\operatorname{div}\left(\frac{Du^\varepsilon}{\sqrt{\varepsilon^2+|Du^\varepsilon|^2}}\right) = -(\varepsilon^2+|Du^\varepsilon|^2)^{-1/(2k)} \quad \text{in } \Omega, \qquad u^\varepsilon=0 \quad \text{on } \partial\Omega,$$9

and the flow is the $Du=0$0-gradient flow of the geodesic length $Du=0$1. The schemes incorporate natural boundary conditions for open curves, and stability can be proved for some of them. The paper applies these discretizations to the Angenent metric to compute rotationally symmetric self-shrinkers for mean curvature flow and to geodesic computations relevant for multi-component phase field models (Garcke et al., 2020).

Boundary-value approximation for surfaces requires additional structure. For mean curvature flow with fixed boundaries, one formulation augments the position evolution by evolution equations for the mean curvature $Du=0$2 and the normal vector $Du=0$3: $Du=0$4 The boundary conditions are

$Du=0$5

together with a conormal derivative condition for $Du=0$6. The semidiscrete analysis introduces a nonlinear Ritz projection to handle the normal boundary constraint and proves optimal $Du=0$7 error estimates for surface position, velocity, mean curvature, and normal vector in discrete spaces with piecewise polynomials of degree $Du=0$8. The numerical experiments corroborate these estimates (Ivaniszyn et al., 25 Apr 2025).

These variational and geometric schemes differ from level-set and diffuse-interface methods in that the primary unknown is the evolving discrete geometry itself. Their main structural assets are discrete energy decay, explicit treatment of boundary conditions, and direct access to tangential mesh motion.

5. Varifolds, point clouds, Brakke-type equalities, and surgery limits

Approximation by varifolds extends mean curvature flow to objects for which classical differential geometry is unavailable. A kernel-based approximate mean curvature for a varifold $Du=0$9 is defined by

$$V_h=\{v\in C^0(\bar\Omega_h): v|_{\partial\Omega_h}=0,\ v|_T \text{ linear}\},$$0

The 2025 generalization classifies linear operators $$V_h=\{v\in C^0(\bar\Omega_h): v|_{\partial\Omega_h}=0,\ v|_T \text{ linear}\},$$1 for which this approximation converges to the classical mean curvature and identifies others for which the limit is zero; it also relaxes regularity assumptions to integral varifolds with unit density and $$V_h=\{v\in C^0(\bar\Omega_h): v|_{\partial\Omega_h}=0,\ v|_T \text{ linear}\},$$2, $$V_h=\{v\in C^0(\bar\Omega_h): v|_{\partial\Omega_h}=0,\ v|_T \text{ linear}\},$$3, provided the singular part of the first variation vanishes. The same framework is extended to an approximate second fundamental form. For point-cloud motions driven by this approximate curvature, the paper proves new comparison principles, including interior and exterior sphere barriers in the continuous setting and barrier statements for implicit and explicit time-discrete schemes, with the sphere radius obeying

$$V_h=\{v\in C^0(\bar\Omega_h): v|_{\partial\Omega_h}=0,\ v|_T \text{ linear}\},$$4

in the continuous comparison model (Sagueni, 8 Sep 2025).

A more global construction defines approximate mean curvature for a general varifold through a mollified first variation,

$$V_h=\{v\in C^0(\bar\Omega_h): v|_{\partial\Omega_h}=0,\ v|_T \text{ linear}\},$$5

and then advances the varifold by iterated push-forwards

$$V_h=\{v\in C^0(\bar\Omega_h): v|_{\partial\Omega_h}=0,\ v|_T \text{ linear}\},$$6

For fixed $$V_h=\{v\in C^0(\bar\Omega_h): v|_{\partial\Omega_h}=0,\ v|_T \text{ linear}\},$$7, the time-discrete flow converges as the time step tends to zero to a unique continuous approximate mean curvature flow $$V_h=\{v\in C^0(\bar\Omega_h): v|_{\partial\Omega_h}=0,\ v|_T \text{ linear}\},$$8. This limit is stable, satisfies a Brakke-type equality, and has monotone mass decay. After coupling with the canonical time measure and letting $$V_h=\{v\in C^0(\bar\Omega_h): v|_{\partial\Omega_h}=0,\ v|_T \text{ linear}\},$$9, one obtains a spacetime limit measure with bounded generalized mean curvature; under an additional rectifiability assumption, that limit is a spacetime Brakke flow (Buet et al., 1 Oct 2025).

Related consistency results have been established for volumetric varifolds. If $u_h^\varepsilon\in V_h$0 is a volumetric discretization of a smooth flow and $u_h^\varepsilon\in V_h$1 its approximate mean curvature, then an integral Brakke approximate equality controls the deviation from the Brakke identity by terms involving the bounded-Lipschitz distance $u_h^\varepsilon\in V_h$2, the test function norms, and the scale-discretization coupling

$u_h^\varepsilon\in V_h$3

This gives a quantitative measure of consistency between volumetric discretization and mean curvature evolution in weak form (Sagueni, 8 Sep 2025).

Approximation by smooth flows with surgery belongs to the same weak-limit landscape, although the mechanism is topological rather than kernel-based. For 2-convex hypersurfaces, mean curvature flow with surgery exists in all $u_h^\varepsilon\in V_h$4; surgeries occur on strong necks, discarded components are geometrically controlled, and as the threshold for discarding high-curvature components tends to infinity the surgical flows converge to the level-set flow (Haslhofer et al., 2014). A later construction shows that weak mean curvature flows with only spherical and neck-pinch singularities can be approximated by smooth flows with finitely many surgeries, with Hausdorff convergence everywhere and smooth convergence away from the singular set (Daniels-Holgate, 2021).

6. Structural distinctions, related flows, and recurrent misconceptions

Approximation theory interacts with a broader question: whether mean curvature flow can itself be realized as a genuine gradient flow on a nondegenerate metric space of hypersurfaces. One related construction is uniformly compressing mean curvature flow, designed to preserve uniform gridpoint density. In the curve case it evolves by

$u_h^\varepsilon\in V_h$5

where $u_h^\varepsilon\in V_h$6 and $u_h^\varepsilon\in V_h$7 is a Lagrange multiplier solving an elliptic equation. The flow is presented as a formal negative gradient flow of the length functional on the manifold of immersed constant-speed curves with fixed center of mass, viewed as a submanifold of Wasserstein space endowed with Otto’s Riemannian structure. In one dimension the paper proves local well-posedness, exponential convergence to a circle near the stationary state, and global weak solutions for arbitrary Lipschitz initial data (Shi et al., 2017).

However, exact mean curvature flow is not a gradient flow on either of two nondegenerate metric structures studied for simple closed plane curves: the uniformness-preserving metric of Shi and Vorotnikov and the curvature-weighted Michor-Mumford metric. The paper establishes this by showing the necessary conservativity conditions fail through explicit counterexamples. The consequence is that several successful “gradient-flow approximations” of mean curvature motion should be interpreted as approximations or alternative projected flows, not as exact metric reformulations of classical mean curvature flow (Huang, 2022).

Another distinction concerns nonlocality. A non-local mean curvature flow based on the first variation of

$u_h^\varepsilon\in V_h$8

is not the classical local flow, although $u_h^\varepsilon\in V_h$9 $$\int_{\Omega_h} \frac{Du_h^\varepsilon\cdot Dv_h}{\sqrt{\varepsilon^2+|Du_h^\varepsilon|^2}}\,dx = \int_{\Omega_h} (\varepsilon^2+|Du_h^\varepsilon|^2)^{-1/(2k)} v_h\,dx \qquad \forall v_h\in V_h.$$0-converges to the perimeter as $$\int_{\Omega_h} \frac{Du_h^\varepsilon\cdot Dv_h}{\sqrt{\varepsilon^2+|Du_h^\varepsilon|^2}}\,dx = \int_{\Omega_h} (\varepsilon^2+|Du_h^\varepsilon|^2)^{-1/(2k)} v_h\,dx \qquad \forall v_h\in V_h.$$1. Its analysis required a viscosity-solution framework based on continuity with respect to Kuratowski convergence of level sets and an exact flow constructed by minimizing movements. Numerical experiments reported qualitative differences from standard mean curvature flow, especially the preservation of small oscillations (Chambolle et al., 2012).

A recurrent misconception is therefore that all approximation schemes are merely computational surrogates for the same geometric object. The literature shows a sharper division. Some methods converge to viscosity solutions of the classical or approximated level-set PDE, some to Brakke flows, some to smooth pre-singular evolutions, and some define alternative geometric motions with improved numerical or metric properties. Analytical estimates on the classical flow remain relevant in this setting: pointwise control of $$\int_{\Omega_h} \frac{Du_h^\varepsilon\cdot Dv_h}{\sqrt{\varepsilon^2+|Du_h^\varepsilon|^2}}\,dx = \int_{\Omega_h} (\varepsilon^2+|Du_h^\varepsilon|^2)^{-1/(2k)} v_h\,dx \qquad \forall v_h\in V_h.$$2 in terms of initial geometry and $$\int_{\Omega_h} \frac{Du_h^\varepsilon\cdot Dv_h}{\sqrt{\varepsilon^2+|Du_h^\varepsilon|^2}}\,dx = \int_{\Omega_h} (\varepsilon^2+|Du_h^\varepsilon|^2)^{-1/(2k)} v_h\,dx \qquad \forall v_h\in V_h.$$3 bounds yields extension theorems and blowup-rate information, which is pertinent when assessing the pre-singular regime of high-accuracy approximations (Wang, 2021).