PaperlayerView exact revision on arXiv

Source-linked AI summary

Variational Derivation of the Weak Form for Viscous Flows Involving Free Surfaces with Implications for Geodynamic Simulations

Karim Norouzi Moghanjoghi, Javier Garcia Pintado, Daniel R. Shapero, Marta Perez-Gussinye

arXiv:2610.12152v1math.NAcs.NAphysics.flu-dynphysics.geo-ph

TL;DR

Free-surface models face time-step restrictions from oscillations, while commonly used stabilization terms have lacked a first-principles derivation. The paper derives the weak form for a deforming domain through variational calculus, finding gravity and viscous boundary contributions. Its benchmarks are consistent with prior stabilization approaches, and the free surface makes the inertia-free problem time-dependent rather than steady-state.

  • Problem

    Larger time steps in free-surface geodynamic models can trigger oscillations and numerical blow-up, while stabilization methods have largely been derived heuristically.

    Source: §1 · p. 2

  • Method

    The authors use an energy-minimization formulation and perturb both the state variables and the deforming domain to derive the weak form and its boundary contributions.

    Source: §1 · p. 3, §4.2 · pp. 8, 11

  • Results

    The formulation yields gravity and viscous boundary terms; benchmark solutions are close to prior stabilization methods, and the free surface makes the inertia-free problem time-dependent rather than steady-state.

    Source: §7 · p. 17, §8 · pp. 19–20

  • Takeaways & Limitations

    The derivation gives a mathematical foundation for gravity stabilization and identifies a viscous boundary contribution that can be significant beyond the simple test cases studied.

    Source: §8 · p. 19

    A comprehensive empirical assessment of the viscous stabilization term is outside the study’s scope, and meaningful validation requires systematic analysis across time-integration schemes.

    Source: §8 · p. 20

Abstract

from arXiv · show

Geodynamic simulations involve solving the Stokes equations, where boundary conditions play a crucial role in governing internal dynamics and surface topography evolution. A common approach is to model the interface in contact with air as a free surface. This interface is stress-free and deforms in response to internal force balances, which is essential for understanding processes ranging from rifted margins to subduction. However, a free surface can introduce instabilities and oscillations in the numerical solution, commonly known as the "drunken sailor" instability, forcing the use of small time steps. Geodynamic models therefore often employ stabilization terms based on the gravitational effect of free-surface motion. To date, these terms have been derived primarily on a heuristic basis. Here, we apply variational calculus to the energy-minimization formulation of incompressible Stokes flow. Using the Gateaux derivative, we derive the complete weak form for a deforming free-surface domain and show that the stabilization terms arise naturally from the free-surface assumption. We consider the linear-viscosity case and identify two boundary-integral contributions: one associated with internal viscous deformation and the other with gravitational body force. Two benchmark cases demonstrate the consistency of the derived formulation with previously proposed stabilization approaches. Our results provide a mathematical foundation for the commonly used gravity stabilization and identify an additional viscous boundary contribution.

1. Introduction

Free-surface boundary conditions affect geodynamic dynamics and topography, but including a free surface restricts time steps because larger steps can trigger oscillations and numerical blow-up. The paper addresses the lack of a first-principles derivation for stabilization by deriving a weak form that includes gravity and viscous terms.

1. Introduction

  • The deforming free surface responds to internal stresses and body forces, and its boundary condition can affect slab-subduction behavior.The cited studies distinguish a true free surface from free-slip conditions and sticky-air approximations.Source: §1 · p. 2
  • Larger time steps can make the velocity field oscillate and cause numerical blow-up, the “drunken sailor” effect.Sticky air, interface capturing, and implicit methods have been explored to mitigate the issue.Source: §1 · p. 2
  • Existing stabilization approaches include gravity-based and other methods, but they have been proposed heuristically rather than derived from first principles.Kaus et al. use a Taylor expansion with a gravity term and a time-step parameter; some implicit approaches add an equation to the system.Source: §1 · pp. 2–3
  • An energy-minimization derivation for a deforming free-surface domain yields a weak form with viscous and gravity stabilization terms, with the gravity term coinciding with Kaus et al..The formulation supplies a broader weak form arising from spatial variation of the domain.Source: §1 · p. 3

2. Governing equations

The section sets out the incompressible Stokes problem and identifies its unknown fields, stress, viscosity, and boundary-condition setting.

2. Governing equations

  • The Stokes formulation is posed for velocity and pressure fields on a bounded domain, with the governing equations introduced for incompressible flow.The passage names velocity u and pressure p as the unknowns.Source: §2 · p. 3
  • The body force is usually gravity, and the Cauchy stress tensor is decomposed into deviatoric and volumetric (pressure) parts.The effective viscosity is positive, and the strain-rate tensor is the symmetric gradient of velocity.Source: §2 · p. 3
  • The strong form uses Dirichlet and nonempty Neumann boundaries; the study assumes zero traction, so no pressure null space is present.Dirichlet boundaries specify velocity, while Neumann boundaries specify traction.Source: §2 · pp. 3–4

3. Derivation of the standard weak form

The standard weak form is built by testing the Stokes equations, integrating over the domain, and applying integration by parts to expose the boundary contribution.

3. Derivation of the standard weak form

  • The weak-form derivation multiplies PDE terms by test functions and integrates them over the domain to form a weighted residual.The passage introduces Lebesgue and Sobolev spaces for the formulation.Source: §3 · p. 4
  • The derivation defines Sobolev spaces for vector-valued functions and admissible velocity test and trial functions.The velocity and pressure test spaces are identified as V0 and Q = L2(Ω).Source: §3 · p. 4
  • Testing the momentum equation and integrating by parts converts the stress-divergence term into a boundary flux contribution.The weak form retains surface traction as a Neumann boundary term and is equivalent to the variational derivation in Boffi et al..Source: §3 · pp. 4–5
  • The weak formulation also tests incompressibility with a pressure test function; tensor colon denotes the Frobenius product, and strain rate is symmetric.These properties support the stress and strain-rate expressions in the weak form.Source: §3 · p. 5

4. Variational approach

The section derives Stokes weak forms by energy variation, first for a fixed domain and then for a deforming free-surface domain. The deforming-domain derivation yields boundary contributions that include gravity and viscosity.

4.1.1. Application to Stokes equation.

  • On a fixed domain, energy minimization and the Gâteaux derivative yield the classical mixed Stokes formulation, with pressure variation enforcing incompressibility.The fixed-domain setting perturbs only state variables, so no geometric terms appear.Source: §4.1 · pp. 6–7, §4.1.1 · p. 7

4.2. Variational approach for deforming domains

  • For a deforming domain, the first variation separates intrinsic changes in state variables from geometric changes induced by domain perturbation.The derivation retains first-order contributions and considers linear viscosity for tractability and comparison with prior approaches.Source: §4 · p. 6, §4.2 · p. 8
  • Domain variations reduce to a boundary integral, while the resulting expression retains the fixed-domain variation when domain perturbation vanishes.The derivation treats domain and state-variable variations separately before applying the divergence theorem.Source: §4.2 · p. 11

4.2.1. Application to Stokes equation with free surface.

  • The free-surface boundary integral contains gravity and viscous contributions; omitting the viscous term when boundary effects are negligible recovers commonly used gravity stabilization.The gravity term coincides with previously proposed stabilization, while the viscous boundary term was not previously considered.Source: §4.2.1 · pp. 12–13
  • As Δt → 0, the boundary term vanishes and the formulation reduces to the standard steady-state Stokes weak form; for larger time steps it accounts for finite free-surface deformation.The authors describe the formulation as consistent at both small and large time steps.Source: §4.2.1 · p. 12

5. Dimensional analysis of the significance of terms

The section nondimensionalizes the functional to compare viscous and gravitational contributions. Its scaling and dimensional analysis provide a first-order estimate of their relative importance in geodynamic systems.

5. Dimensional analysis of the significance of terms

  • The derivation factors the original functional into a dimensional prefactor and a purely dimensionless functional.The section introduces dimensionless variables, then defines I* before expressing the original functional in that form.Source: §5 · pp. 13–14
  • The functional scales with ρ, U, and L, while I* is dimensionless, allowing estimates of the relative viscous and gravitational contributions in geodynamic problems.The section uses dimensionless variables and parameters to assess the terms in the weak-form boundary integral.Source: §5 · pp. 13–14
  • Re/Fr2 = ρgL2/(µU) quantifies gravitational forcing relative to viscous resistance in the Stokes regime.Table 1 presents characteristic scales for representative geodynamic systems, with final-column values rounded to the nearest order of magnitude.Source: Table 1 · p. 14
  • Dimensional analysis indicates that the viscous boundary term matters in rifting and plate subduction, and should not be neglected in mantle convection or ice-sheet and glacial systems.For lower-viscosity lava flow and magma ascent, the viscous term can be ignored; the analysis is a preliminary estimate, not a definitive criterion for model evolution.Source: §5 · p. 14

6. Numerical Test Cases

Two test cases, ordered by increasing physical and numerical complexity, assess the accuracy, stability, and applicability of the free-surface stabilization formulations. They examine free-surface behavior in an unstable layering and in a rising-plume setting.

6.1. Rayleigh–Taylor Instability

  • A denser, more viscous layer over a less dense, less viscous layer forms the gravitationally unstable Rayleigh–Taylor benchmark.It is a standard test of free-surface method stability and accuracy.Source: §6.1 · p. 15
  • The Rayleigh–Taylor test examines whether stabilization suppresses free-surface oscillations while preserving physical instability growth across viscosity contrasts and time-step sizes.Large time steps may produce oscillatory surface motion and eventual divergence from the physically evolving solution.Source: §6.1 · p. 15

6.2. Rising Plume

  • The rising-plume benchmark uses a true free surface tracked in a Lagrangian framework, assessing surface deformation without an artificial low-viscosity layer.The setup follows the idealized plume benchmark used by Crameri et al. to investigate surface-topography errors associated with the sticky-air approximation.Source: §6.2 · p. 15
  • The plume test compares predicted surface evolution and topography while varying lid viscosity; constant viscosity within each material isolates stabilization effects from nonlinear rheology.The buoyant plume interacts with a deformable surface, whose topography depends on the relative viscosities of the plume and surrounding layers.Source: §6.2 · p. 16, Figure 3 · p. 16

7. Results

The variational formulation’s full stabilization is consistent with gravity-only stabilization and closely tracks the small-timestep reference in both benchmarks. Across the Rayleigh–Taylor tests, the methods remain consistent across lid viscosities, while the rising-plume results show no significant topography errors from either method.

7. Results

  • Both stabilization methods closely track the dt = 100-year reference in the Rayleigh–Taylor and rising-plume benchmarks; the full and gravity-only results are nearly identical.The Rayleigh–Taylor comparison uses dt = 5000 years, and the rising-plume comparison also reports dt = 5000 years.Source: §7 · pp. 17, 19, Figure 4 · p. 18, Figure 6 · p. 19
  • Across tested time steps and lid viscosities, both methods perform consistently; Rayleigh–Taylor MAE peaks at lower viscosity, where surface deformation is greatest.The larger discrepancies at low viscosity remain below typical error margins, and observed variance is lower than reported for sticky-air methods.Source: §7 · p. 17, Figure 5 · p. 18
  • As lid viscosity rises, maximum elevations agree better and minimum elevations vary less; total relief peaks at µ = 1023 Pa·s and is comparable at µ = 1022 and 1024 Pa·s.Neither stabilization method introduces significant errors into predicted topography.Source: §7 · p. 19
  • “Full stabilization” combines viscous and gravity effects, whereas “Gravity-only stabilization” applies only the gravity term in the boundary integral.These labels distinguish the proposed formulation from the gravity-only comparison in both test cases.Source: §7 · p. 19

8. Discussion and Conclusion

The variational formulation produces gravity and viscous boundary terms, with their relative importance tied to Froude and Reynolds numbers. Its scope is linear-viscosity problems, and systematic empirical assessment of the viscous term remains future work.

8. Discussion and Conclusion

  • The free-surface boundary integral contains gravity and viscous contributions, whose importance scales with the Froude and Reynolds numbers, respectively.The viscous contribution may be significant in many problems, although it is smaller in the simple benchmark cases studied.Source: §8 · p. 19
  • Table 2 reports consistent relief statistics for full and gravity-only stabilization across the tested viscosities in the rising-plume case.The comparison is after 3 Myr with a 5,000-year time step.Source: Table 2 · p. 20
  • The formulation is derived under a linear-viscosity assumption.Source: §8 · p. 20
  • A comprehensive empirical assessment of the viscous stabilization term is outside this study’s scope, and its importance requires dedicated validation.Free-surface dynamics may be sensitive to the time-integration scheme used for the boundary integral; limited testing could be misleading.Source: §8 · p. 20
  • The authors propose controlled comparisons across parameter regimes at identical physical timestamps for unstabilized, gravity-only, and full stabilization.They identify comparisons across multiple settings as necessary to avoid misleading conclusions.Source: §8 · p. 20
  • The free-surface boundary differs from a Neumann boundary because it makes the problem time-dependent even when inertial effects are ignored.The authors suggest this foundation may inform time-stepping methods for models involving surface and subsurface interaction.Source: §8 · pp. 20–21

Appendix A. Supplementary Material

The supplementary material provides numerical implementation details and quantitative verification, complementing the main text’s governing equations and physical model.

Appendix A. Supplementary Material

  • The supplementary material adds details on the numerical workflow, ALE mesh movement, remeshing, and verification benchmarks.Source: Appendix A · p. 22

Appendix A.1. Numerical Workflow and ALE Time Integration

The ALE time-stepping procedure separates physical flow from mesh motion and updates the computational mesh after solving the Stokes problem.

Appendix A.1. Numerical Workflow and ALE Time Integration

  • At each time step, the Stokes problem is solved on the current mesh, and its physical velocity then determines selected mesh-boundary motion.The physical velocity evolves modeled material, while mesh velocity moves the computational grid and helps retain element quality.Source: Appendix A.1 · p. 22
  • The computational state at time tn includes mesh coordinates, physical fields, and solver and function-space definitions.The time-stepping procedure is summarized in Algorithm 1.Source: Appendix A.1 · p. 22
  • The mesh coordinates are updated explicitly, mesh quality is checked, and a mesh is regenerated when its minimum quality measure falls below a threshold.The workflow then advances the time controller; reconstruction of mesh or function spaces triggers their rebuilding and solver reinitialization.Source: Appendix A.1 · p. 22

Appendix A.2. Mesh Movement

Mesh motion is computed separately from the Stokes flow through an elliptic smoothing problem, with prescribed boundary movement extended into the domain.

Appendix A.2. Mesh Movement

  • The mesh velocity is obtained by solving a separate linear elliptic smoothing problem with conditions for moving material-following boundaries and fixed boundaries.Source: Appendix A.2 · p. 23
  • Harmonic extension propagates prescribed boundary motion into a smooth interior displacement field, after which mesh coordinates advance explicitly.Source: Appendix A.2 · p. 23
  • Separating the mesh-movement solver from the Stokes solver allows changes to mesh motion or boundary tracking without changing the physical flow formulation.Source: Appendix A.2 · p. 23

Appendix A.3. Mesh Quality and Remeshing

The appendix describes mesh-quality monitoring and remeshing, then verifies the implementation with complementary smooth and heterogeneous analytical benchmarks. Together, the tests provide a progressive framework for assessing accuracy and numerical reliability.

Appendix A.3. Mesh Quality and Remeshing

  • The global mesh-quality check uses minimum element angle, and excessive accumulated distortion triggers remeshing while retaining the represented geometry and interfaces.After remeshing, the finite-element spaces and solution fields are rebuilt and the numerical solution is transferred to the new mesh.Source: Appendix A.3 · p. 23

Appendix A.5. Manufactured-Solution Benchmark

Appendix A.6. SolVi Viscous-Inclusion Benchmark

  • In SolVi, the numerical solution reproduces the principal pressure extrema and symmetry, while the domainwide fields show perturbations from the high-viscosity inclusion.Figures A.9 and A.10 show the centerline pressure comparison and domainwide pressure and velocity fields.Source: Appendix A.6 · p. 27, Figure A.9 · p. 27, Figure A.10 · p. 27
  • In the SolVi benchmark, the largest discrepancies occur at the circular viscosity interface, where the discontinuous material coefficient reduces solution regularity.The benchmark tests the formulation against a strong viscosity contrast and sharp internal interface.Source: Appendix A.6 · pp. 26–27

Appendix A.4. Numerical Verification and Benchmarking

  • The verification procedure pairs a smooth known solution for baseline accuracy with SolVi’s strong viscosity contrast and discontinuous interface as a more demanding test.Pressure errors are evaluated after both numerical and analytical pressures are normalized to zero mean, removing the arbitrary pressure gauge.Source: Appendix A.4 · p. 24, Appendix A.7 · p. 28

Appendix A.7. Convergence and Verification Procedure

  • The two benchmarks, combined with ALE mesh movement and remeshing, form a reproducible framework for assessing simulation accuracy and numerical reliability.The manufactured solution establishes baseline accuracy; SolVi evaluates performance with strong material heterogeneity and a sharp interface.Source: Appendix A.7 · p. 28
Loading 2610.12152v1…