Benchmarking Stokes Flow in Annulus and Spherical-Shell Geometries

Photo by the blowup / Unsplash

This post introduces a suite of annulus and spherical-shell Stokes benchmarks reproduced using Underworld3 1. These benchmark problems have previously been implemented in several numerical codes, but here we bring them together within a single framework to highlight both the strengths and practical challenges of curved-domain finite-element modelling. The post is aimed at researchers working on geodynamics numerical modelling. Rather than focusing on heavy mathematical derivations or low-level implementation details, the goal is to provide an intuitive and practical guide to what each benchmark is designed to test, what numerical behaviour should be expected, and how to interpret the results. By the end, readers should have a clear understanding of the purpose of these benchmarks and the key ideas behind verifying Stokes flow in curved geometries.

What Is Benchmarking and Why Is It Important?¶

Benchmarking is the process of testing a numerical method or software implementation against problems with known analytical solutions or well-established reference results. In computational geodynamics, benchmarks are essential because they help verify that a code correctly solves the governing equations before it is applied to complex Earth-science problems where the true solution is unknown. A good benchmark does more than produce a visually reasonable result; it tests numerical convergence behaviour, boundary-condition implementation, mesh geometry, and solver robustness under controlled conditions. Curved-domain benchmarks such as annulus and spherical-shell Stokes problems are particularly important because they expose numerical challenges that do not appear in simple Cartesian geometries. Successfully reproducing benchmark results therefore builds confidence that a numerical framework can model real-world problems with reliable and physically meaningful approximations.

What Are the Stokes Equations?¶

The Stokes equations describe the slow, viscous flow of fluids in situations where inertial forces are negligible compared with viscous forces. In geodynamics, this approximation is widely used because rocks in the Earth’s mantle deform extremely slowly over geological timescales and behave like highly viscous fluids. The incompressible Stokes equations are written as

Here $\mathbf{u}$ is velocity, $p$ is pressure, $\eta$ is viscosity, $\rho$ is density, and $\mathbf{g}$ is gravity. The first equation represents conservation of momentum, balancing viscous stresses, pressure gradients, and body forces. The second equation enforces mass conservation through incompressibility. Solving these equations allows us to model mantle convection, lithospheric deformation, subduction, and many other large-scale Earth processes. Although the equations appear compact, solving them accurately in curved geometries with complex boundary conditions and variable material properties is computationally challenging, which is why benchmark problems are important.

Why Do We Need Curved Geometries?¶

Curved geometries are important in geodynamics because the Earth itself is curved. Many large-scale Earth processes, such as mantle convection, subduction, plume dynamics, and lithospheric deformation, occur within spherical or shell-like domains rather than simple rectangular boxes. While Cartesian geometries are useful for developing intuition and testing numerical methods, they cannot fully represent radial gravity, curved boundaries, or global-scale flow patterns. Curved-domain models also introduce additional numerical challenges, including geometric approximation errors, coordinate transformations, and the accurate implementation of free-slip or zero-slip boundary conditions on non-planar surfaces. As a result, annulus and spherical-shell benchmarks provide a more realistic and demanding test of Stokes solvers.

Benchmark Suite¶

In this work, we reproduce four widely used Stokes benchmark suites in curved geometries using Underworld3 1. The first is the Thieulot–Puckett annulus benchmark, which provides a smooth analytical solution in an annulus and is mainly used to test optimal finite-element convergence behaviour 2. The second is the Kramer annulus benchmark, which extends the problem to include both smooth volumetric forcing and singular delta-function forcing on an internal interface, together with free-slip and zero-slip boundary conditions 3. We then consider the spherical-shell counterparts of these problems: the Thieulot spherical benchmark, which tests smooth Stokes flow in spherical geometry with both constant and radially varying viscosity 4, and the Kramer spherical benchmark, which again introduces internal interface forcing and reduced solution regularity 3. Together, these four benchmark suites test curved geometries, pressure treatment, mesh approximation, boundary-condition implementation, smooth and singular forcing, and convergence behaviour in both two- and three-dimensional Stokes flow.

What Do We Measure? Error Quantification¶

The benchmark comparisons use $L_2$-norm errors because they give a single quantitative measure of the difference between the numerical solution and the analytical solution over the whole domain. For a computed field $q_h$ and analytical field $q^{*}$, the absolute volume error is

The corresponding relative volume error is

The relative form is useful when comparing velocity and pressure errors across different benchmark cases because it normalises the error by the size of the analytical solution. For pressure, the numerical and analytical fields are first compared in the same pressure gauge, since incompressible Stokes pressure is determined only up to an additive constant.

Boundary pressure errors are measured separately on the inner and outer surfaces. For a boundary $\Gamma \in \{\Gamma_{\mathrm{inner}}, \Gamma_{\mathrm{outer}}\}$, the absolute pressure-trace error is

Where the analytical boundary pressure has a nonzero $L_2$ norm, the relative boundary pressure error is

Convergence is measured by comparing errors across successively refined meshes:

Here $h$ is the characteristic cell size. If the mesh is uniformly refined so that $h_2 = h_1/2$, this becomes

The expected behaviour is simple: as $h$ decreases, the error should decrease. On a log-log convergence plot, a method with error proportional to $C h^{r}$ appears approximately as a straight line with slope $r$. For smooth Stokes solutions, stable mixed finite-element pairs have well-defined optimal convergence expectations; for example, Taylor–Hood $P_2 \times P_1$ commonly gives third-order velocity and second-order pressure convergence in the volume $L_2$ norm when the geometry and solution are sufficiently smooth 5. Singular forcing cases are different: the solution is less regular near the internal interface, so reduced convergence rates are expected even when the solver is implemented correctly.

What Do We Observe? Results of Benchmarks in Underworld3¶

Volumetric Convergence¶

The benchmark results show that Underworld3 reproduces the expected convergence behaviour for both annulus and spherical-shell Stokes problems. For smooth analytical solutions, such as the Thieulot annulus and spherical benchmarks, the Taylor–Hood $P_2 \times P_1$ discretisation achieves close to the theoretically expected convergence rates, with approximately third-order velocity convergence and second-order pressure convergence in the volume $L_2$ norm. Higher-order element pairs further improve accuracy, while lower-order discretisations show the expected reduction in convergence order.

For the Kramer benchmarks with smooth forcing, the velocity convergence is closer to second order because the curved geometry is represented using linear meshes, consistent with the original benchmark studies 3. In the delta-function forcing cases, the convergence rates reduce significantly because the internal singular interface lowers the regularity of the analytical solution. Overall, the volumetric results demonstrate that Underworld3 captures both optimal convergence behaviour for smooth problems and the expected degradation in accuracy for singular forcing cases.

Boundary Pressure Convergence¶

Boundary diagnostics provide a more local and often more sensitive measure of solver accuracy in curved geometries. The volume pressure norm measures an error over the full domain, but boundary pressure errors isolate the pressure trace on the surfaces where boundary conditions, radial normal stresses, and traction-related diagnostics are evaluated. This is useful because geometric approximation errors and boundary quadrature errors can be more visible on curved inner and outer boundaries than in a volume-averaged metric.

For the Thieulot–Puckett annulus benchmark, the $P_2 \times P_1$ boundary pressure errors decrease close to the expected $O(h^2)$ trend. The inner and outer boundaries do not have identical error constants, which is acceptable: the two curves have different radii, different geometric representation errors, and different boundary integration paths.

For the Kramer annulus benchmark, the boundary pressure errors converge close to second order for the smooth cases. The delta-function cases also show stronger boundary pressure convergence than the corresponding volume pressure norm because the reduced regularity is localized at the internal interface rather than on the annulus boundaries.

The Thieulot spherical-shell benchmark gives the corresponding smooth three-dimensional boundary test. The boundary-pressure panel shows that the absolute pressure-trace errors for both $m = -1$ and $m = 3$ decrease systematically with refinement. These errors are shown together with radial normal-stress errors because the stress diagnostic combines pressure and velocity-gradient errors at the same curved spherical boundaries.

For the Kramer spherical-shell benchmark, the current boundary figure reports boundary velocity and radial normal-stress convergence for the free-slip delta-function case. A separate boundary pressure-trace convergence plot is not part of the current spherical Kramer article outputs. The radial normal-stress diagnostic is included here because it is pressure-sensitive: $\sigma_{rr}$ contains the pressure contribution and therefore tests pressure recovery together with velocity-gradient recovery on the spherical boundaries.

Overall, the boundary diagnostics complement the volume $L_2$ pressure errors. Smooth benchmark cases recover approximately second-order pressure-trace convergence for the $P_2 \times P_1$ pair, while singular forcing primarily affects the volume pressure norm through reduced regularity at the internal interface. Persistent differences between inner- and outer-boundary error constants should therefore be interpreted as geometric and boundary-evaluation effects rather than as a failure of pressure normalisation.

The full details are in the article PDFs, which can be viewed directly from the benchmark figures and articles directory on GitHub:

• Thieulot–Puckett annulus benchmark

• Kramer annulus benchmark

• Thieulot spherical-shell benchmark

• Kramer spherical-shell benchmark

Those articles contain the analytical fields, discretisation choices, convergence tables, boundary diagnostics, and full reference details.