Volume of Fluid method applied to free surface boundaries in numerical geodynamic models

Abstract. Tracking the evolution of interfaces in numerical geodynamic models, particularly the rock-air/water boundary when imposing free surface boundary conditions, presents significant computational challenges. Traditional marker-in-cell methods, while widely used in mantle convection codes like StagYY, suffer from inherent noise due to random tracer distributions and require high tracer densities for accurate topography resolution, leading to substantial computational costs. This study presents the implementation of the Volume of Fluid (VOF) method into StagYY using the open-source gVOF package, offering a volume-conservative alternative for interface tracking in geodynamic simulations. Two implementations are developed and tested: a basic method with one VOF cell per computational cell, and an improved method utilizing subcells with full coupling to the Stokes solver through density and viscosity fields. Both methods employs the CLCIR (Conservative Level Contour-based Interface Reconstruction) reconstruction scheme combined with FMFPA (Face-Matched Flux Polyhedron Advection) for optimal accuracy. A consolidation approach for translating VOF fields to Eulerian surface locations ensures volume conservation and numerical stability. Benchmarking in 2D and 3D across multiple geometries demonstrates that while the basic VOF implementation provides adequate results for simple 2D cases, the improved method is essential for accurate tracking in complex flows and three-dimensional applications. The improved VOF method successfully eliminates topographic noise inherent in tracer-based approaches and decouples surface quality from tracer density, enabling high-quality results with reduced computational overhead from tracer advection. Despite increased computational costs and memory consumption associated with the method, particularly with the improved implementation, the VOF method offers distinct advantages including explicit volume conservation, sharp interface representation, and versatility for tracking various geophysical boundaries beyond free surfaces. This work establishes a foundation for future applications in sea level modelling, continental margin tracking, and coupled planetary system simulations, advancing toward global-scale biogeodynamic modelling.