See the attached pdf.

This is an interesting and timely paper presenting the new free surface treatment algorithm with the broadly used ALE FEM PIC code Underworld. The new method shows good performance for several test cases compared to two other methods tested with the same code. The paper is suitable for publication after addressing some comments, which are listed below.

Taras Gerya, Zurich, 09.04.2026

Specific comments

Abstract. “The commonly used ’Sticky Air’ method, while practical, suffers from several limitations, including increased computational cost and marker fluctuation issues.” Free surface numerical fluctuation (“drunken sailor instability”) is also characteristic for true Lagrangian free surface treatment (Kaus et al., 2010).

Line 45. “Importantly, the "sticky air" layer is not intended to represent a physical reality; it possesses the same density as air but a viscosity that is on the order

of 10^14 times greater.” Air viscosity is on the order of 10^-5 Pa*s. Commonly used viscosity of sticky air is on the order of 10^17-10^19 Pa*s. Hence, the difference is rather 10^22-10^24 times.

Line 105. “The density of this layer is set to zero, ensuring it exerts no pressure on the actual free surface (the interface between the air and lithosphere).” This is not always the case. Sticky water density is set to 1000 kg/m^3 (e.g., Gerya and Yuen, 2003) to approximate a water-loaded free surface.

Line 125

“(2) Free Surface Resampling

In accordance with the ALE scheme, the x-coordinates Xx in 2D or the x, y-coordinates Xx,y in 3D of the mesh nodes

remain constant. Consequently, we need to resample the vertical coordinates Xz at these specified locations.

(3) Mesh Regridding

130 To achieve a uniform distribution of displacements Dz in the vertical mesh coordinates, we solve Laplace’s equation:

∇2Dz = 0 (7)

The boundary conditions applied here are Dirichlet constraints, which define the top and bottom boundaries as zero and the

internal boundary as new displacement (Dz = Xn+1z −Xnz on Γfs).

Next, we update the vertical mesh coordinates forward in time using displacements determined by the forward Euler scheme:

Xn+1z = Xnz +Dz”

This is confusing. If horizontal coordinates of the topography points do not change then advection of the Lagrangian moving free surface relative to the immobile Eulerian topography points needs to be taken into account when computing new topography elevation Xn+1z. Is this the case with the applied approach?

Line 235. “…and is limited within nine orders of magnitude by applying upper and lower bounds: ηmax = 105η0 and ηmin = 10−4η0, where η0 is a reference viscosity.” Please specify how η0 is defined (from the dislocation creep parameters? Empirically?).

Line 265. “When the free surface is not explicitly tracked using additional tracers, the surface becomes unidentifiable, as shown in Fig. 1e. In such cases, the surface must be tracked via particles representing the top of the solid or the interface between rock

and air, or through an averaged interface based on volume ratios (Deng et al., 2024). Using extra particles to trace the surface, common in this study, often results in a rough interface with undesired spatial fluctuations as discussed in Crameri et al. (2012). These fluctuations arise because the distance between markers and the interface is finite and irregular, leading to small velocity variations during advection. Such fluctuations can be mitigated by employing finer vertical spacing in the computational mesh or by utilizing marker

chains or level-set methods to more accurately assign viscosity and density to nodal points.” In case of the Eulerian MIC approach, somewhat better topography surface representation could be also achieved by: (i) interpolating material type (0=Air, 1=rocks) from Lagrangian markers to Eulerian nodes with distance dependent average followed by (ii) computing position of the topography surface (defined as the isosurface of 0.5) from the Eulearian nodal material type values.

Line 285. “Both the ALE and ALE-IB schemes exhibit excellent agreement in tracking the evolution of the interfaces and the free surface. In contrast, the Eulerian scheme displays significant fluctuations in both the free surface and the lithosphere/asthenosphere interface, along with asymmetric features, especially in the interface’s depth profile (Fig. 6d). The fluctuations observed in the Eulerian approach are likely attributable to the inherent numerical diffusion and irregularities associated with fixed-grid advection, which can cause the interface to oscillate and distort over time.” The poor performance of the Eulerian MIC scheme is somewhat surprising and could be related to the inaccuracy of the applied marker advection scheme (4th order Runge-Kutta, continuity-based advection scheme should typically work better) and/or to the way or recovering the free surface position from markers (Eulearian material type isosurface interpolated from markers should typically work better, see above). Is this specific to the used FEM approach? Does a standard staggered Eulearian grid MIC produce similarly poor results (e.g., freely accessible 2D code example Sticky_air.m from Chapter 8 of Gerya et al., 2019, Second edition)?

https://www.cambridge.org/ch/universitypress/subjects/earth-and-environmental-science/structural-geology-tectonics-and-geodynamics/introduction-numerical-geodynamic-modelling-2nd-edition?format=HB#resources

Line 315. “Our ALE-IB scheme results are more comparable to the free-surface case in the Eulerian scheme reported in Crameri et al. (2017), where a shape-function averaging method was employed in their modelling on all the uppermost rock tracers and the lowermost air tracers. This approach, combined with their sticky air method, yields more accurate surface representations.” Would be good to show their results for comparison on your plots. Will your Eulerian results become more accurate with the use of a similar shape-function averaging method as Crameri et al. (2017) used?

Line 330 “(3) Mesh elements type: In our experiments, all models utilized meshes with Q1dQ0 elements. However, as demonstrated in Thieulot and Bangerth (2022), Q1dQ0 elements tend to be unstable and inaccurate in practice. Consequently, we believe that

higher-order elements such as Q2Q1 and Q2P−1 offer more robust and reliable options for geodynamic simulations, despite their increased implementation complexity and higher computational costs associated with solving the resulting linear systems.” Does Underwold have these higher order elements? If so, presenting some additional tests would be very useful.

Line 350. “Overall, our findings highlight that the ALE-IB scheme not only matches the accuracy of existing methods but also offers significant advantages in stability, robustness, and physical realism. Consequently, it presents a promising alternative to conventional ALE and "sticky air" techniques in the Eulerian scheme, particularly for multi-material near free surface systems and surface process modelling, where precise and stable free surface representation is crucial. This framework paves the way for more reliable and versatile geodynamic simulations, advancing our understanding of Earth’s lithospheric and mantle dynamics with a true free surface.” Advantages over the Eulerian scheme are only demonstrated for the specific scheme explored in Underworld and explored in this study. It performs notably worse than the Eulerian staggered grid MIC scheme of Cramery et. (2017) (see comments to line 315 above).

References

Gerya, T. V. and Yuen, D. A.: Rayleigh–Taylor instabilities from hydration and melting propel ‘cold plumes’ at subduction zones, Earth and Planetary Science Letters, 212, 47–62, https://doi.org/10.1016/S0012-821X(03)00265-6, 2003.

Kaus, B. J. P., Mühlhaus, H., May, D. A. (2010) A stabilization algorithm for geodynamic

numerical simulations with a free surface. Physics of the Earth and Planetary Interiors,

181, 12–20.