DynEarthSol v2.0: An Efficient Explicit Lagrangian Solver for Geodynamics, Surface Processes, and Earthquake-Cycle Dynamics

Shyu, Chase J.; Lee, Sungho; Ding, Xuesong; Keum, Jaeyoon; Tan, Eh; Lavier, Luc L.; Choi, Eunseo

doi:10.5194/egusphere-2026-2922

Preprints

https://doi.org/10.5194/egusphere-2026-2922

Preprints

12 Jun 2026

| 12 Jun 2026

Status: this preprint is open for discussion and under review for Geoscientific Model Development (GMD).

DynEarthSol v2.0: An Efficient Explicit Lagrangian Solver for Geodynamics, Surface Processes, and Earthquake-Cycle Dynamics

Chase J. Shyu, Sungho Lee, Xuesong Ding, Jaeyoon Keum, Eh Tan, Luc L. Lavier, and Eunseo Choi

Abstract. DynEarthSol version 2.0 (DES v2.0) provides a scalable two- and three-dimensional geodynamic modeling platform designed for multiphysical simulations across temporal and spatial scales. This release achieves high-performance execution by leveraging heterogeneous CPU-GPU architectures; and maintains computational consistency on unstructured meshes through a deterministic two-stage update algorithm, eliminating race conditions and ensuring bitwise reproducibility. The updated framework further expands physical realism through a staggered two-way coupling with the goSPL landscape evolution model, which we validated by simulating the development of a low-angle normal fault and the associated geomorphic signatures. Additionally, we integrated a Dieterich-Ruina rate-and-state friction formulation into the existing Mohr-Coulomb elastoplastic structure. By resolving state-variable evolution within finite-thickness shear zones, this approach enables the simulation of complex fault behaviors, from stable creep to stick-slip instabilities, with finite-thickness faults. The new friction model was verified against analytical benchmarks and applied to earthquake-cycle applications. DES v2.0 serves as a robust, GPU-accelerated, and unified platform for investigating the coupled system of tectonic deformation, surface processes, and seismic-cycle dynamics.

Received: 21 May 2026 – Discussion started: 12 Jun 2026

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.

Download & links

Chase J. Shyu, Sungho Lee, Xuesong Ding, Jaeyoon Keum, Eh Tan, Luc L. Lavier, and Eunseo Choi

Status: open (until 24 Aug 2026)

Post a comment Subscribe to comment alert

RC1: 'Comment on egusphere-2026-2922', Anonymous Referee #1, 19 Jul 2026 reply

Shyu and co-authors present implementations of new algorithms in the DynEartSol v2.0 software, which are relevant for the modeling of coupled geodynamic and surface processes in crustal-scale models. In particular they discuss leveraging heterogeneous CPU-GPU architectures, coupling to surface process modeling software, and a new implementation of rate-and-state friction models.

I generally find their manuscript to be well written and the algorithms to be relevant to the field and worth publishing. I have left a number of comments below that I would like the authors to address in order to improve the clarity of the manuscript and the impact of their results. If the authors address these comments I am happy to recommend this manuscript for publication.

Comments:

- line 16: It would be worth mentioning FD and FV approaches here as well. In general the "memory wall" is more related to the existence of an assembled matrix and sparse-matrix vector multiplications (or Krylov solvers with long recurrence relations), not which particular method is used to discretize the equations.

- line 19: this choice of references is not sufficient to claim that "the field" is shifting to matrix-free, explicit formulations, since they only reference the evolution of FLAC. If you want to insist on the "field" and explicit formulations you should include references to other matrix-free explicit methods, such as Lattice Boltzmann methods (e.g. Mora et al 2020, or other works by the authors or others) and other groups making use of similar pseudo-transient methods as FLAC (de Montserrat et al 2026 or other works). However, more broadly even for implicit methods the push to matrix-free methods has the same justification and momentum (e.g. May et al 2015, Clevenger & Heister, 2021, Ilangovan et al 2026). So generally speaking it is more correct to say the field is pushing towards matrix-free methods, independently of whether they approximate matrix-vector products without assembling and storing the matrix, or whether the method completely avoids the use of matrices.

- line 80: Here it would be good to at least mention the stability conditions of the chosen element as well. Linear elements may be a sufficient balance between accuracy and computational efficiency, but depending on their stability they may still not be a good choice for a particular problem. Depending on the exact choice of equations and problem rheology this may or may not be a problem in practice, but it should be discussed in a section about the finite element choice. Also the section could be a bit more clear on which variables are the independent solution variables (I assume force, velocity, and temperature) and which ones are derived (stress, strain, material properties). Also you mention that you use linear basis functions, but from the description I assume they apply to the solution variables, and if stress, strain and material properties are computed at element centroids they are assumed to be cell-wise constant properties, is this correct? If so, please state this in the text. If not, clarify that linear functions are used for all properties.

- Section 3.2, line 133: This sounds as if the GPU part of the code is Nvidia-hardware specific (or was the Nvidia compiler simply used for the specific systems used for testing?). Please clarify in the text if the code for the accelerator is hardware-agnostic and would run on AMD, Intel or other accelerator hardware as well.

- line 189-190: If "profiling confirms" this is the case it would be useful to quantify the amount of wall time spend on CPU vs GPU parts of the algorithm. If this is already computed in the benchmark section later it would be worth referencing the relevant sections here.

- If I/O becomes a bottleneck is there potential for overlap of CPU and GPU operations by putting I/O into a background thread and let is complete on its own terms, while the main computation loop continues?

- line 203: For the benchmarked GPU systems, please also specify in the manuscript which CPUs where used on those systems, as it affects the hybrid CPU/GPU parts of your algorithm. From descriptions in the text below it seems like you used nodes that combined Grace CPU with H100 GPUs, and Epyc 7763 CPUs with A100 GPUs, is this correct? This may affect some of the A100 vs H100 performance comparison (the ones with hybrid algorithms), as the Epyc 7763 is slower than the Grace CPU.

- line 205: "identical configuration, except for mesh and domain sizes" is an unfortunate formulation in this context, because of course if the mesh and domain sizes change, the comparison is no longer fair. I think what you mean is that you performed the identical scaling tests (with variable mesh and domain sizes) on different hardware. Please clarify the text.

- line 217: Grace (Grace) probably should be Grace (CPU)? Also the relevant information for the power law is the exponent of the power law (is it linear, less than linear, something else?), please add it to the text.

- line 218: outperformed the Grace -> outperformed the Grace CPU

- line 229: Again the important information is the (approximate) exponent of the power law and to which model sizes the performance scales, please state it in the text

- line 227: That the EPYC 7763 underperformed the Grace CPU is not surprising given their theoretical performance comparison, memory bandwidth and their relative age (2021 vs 2023). Maybe that is worth mentioning to avoid the impression that AMD Epycs generally underperform the Grace CPU (e.g. the Epyc 9965 from 2024/2025 is reported to have several times the performance of the 7763, and significantly more than the Grace CPU, so the comparison would likely look different). At least specify the CPU model here again.

- Fig. 3: Please double-check that panels b and e are not just accidental copies from panels a and d. They look identical to me, but should not be identical, as a and d should be a sum of all computational steps.

- Section 3.4 in general: One of the main challenges of GPU acceleration for HPC modeling software is the parallellization across distributed memory parallel systems. This is not mentioned here and all tests are performed on one server node, so I assume that DES is purely parallelized for shared memory parallellization? That would be worth pointing out, as it is an important consideration for potential users.

- line 344: This is interesting, it sounds like DES time steps are generally equal or smaller than GoSPL time steps. In other cases of coupling (e.g. ASPECT with Fastscape or Badlands with Underworld) generally the tectonic time step is (often significantly) larger than the surface evolution time step. Is this due to the explicit DES time step? Would the code work correctly if the DES time step is ever significantly larger than the GoSPL time step or is the DES time step to be limited to at most one GoSPL time step?

- line 361: How is the time-averaged deltaX over the coupling period computed? Is there a back up of the node position at the beginning of every coupling period stored? If so, who owns this backup (DES or goSPL)? And how does this work if the DES grid is remeshed during the coupling period? This might be worth adding to the description of this step.

- line 380: You write a single IDW pass is necessary (from goSPL to DES), but step 1 describes that velocity is (optionally) transferred from DES to goSPL via IDW, which would make that 2 IDW passes. This seems inconsistent. Is velocity transfer usually not required and therefore in a typical model only one IDW pass is necessary?

- line 382: It would be good to mention if the size of the padded area sets an upper limit on the maximal extensional strain a DES model can accumulate during the model runtime. Or can the goSPL domain be reset after a certain amount of extension has accumulated?

- line 418: "between" -> "compared to" the DES-only case?

- line 455: This section needs a sentence summarizing the result of the comparison. Is the comparison between DES-goSPL and ASPECT-FASTSCAPE satisfactory? Are the remaining differences likely caused by small differences in model setup (e.g. different bottom boundary condition and/or minor implementation differences)? Or may they be caused by some fundamental disagreement? I would be particularly curious about the authors opinion on the different roughness of the topography between the two models. Are these models run on similar resolutions? Is the higher roughness of the DES-goSPL model caused by the less structured mesh of DES, or does ASPECT-Fastscape seem to be more diffusive? I understand that a conclusive answer to these questions may be beyond the scope of the paper.

- line 592: compared with analytical solution -> compared with an analytical solution

- line 622: I appreciate this benchmark case with analytical solution, which forms a great basis for intercode comparisons! That being said, "are visually indistinguishable" can only be the very first step in a comparison between numerical and analytical benchmark solution. Since you have both solutions available, could you (a) at a minimum list the typical (or final) quantitative relative error between numerical and analytical solution, and (b) consider if it is worth plotting the relative error for the different cases in Fig. 10 in addition to the absolute stress values? That would make it easier to use this benchmark in the future to compare accuracy across different software implementations.

- line 674: "can handle the earthquake with ..." either "an earthquake ..." or "can handle earthquakes ..."

- line 675; "elastic wave" -> "elastic waves"

References:

Peter Mora, Gabriele Morra, David A Yuen, A concise python implementation of the lattice Boltzmann method on HPC for geo-fluid flow, Geophysical Journal International, Volume 220, Issue 1, January 2020, Pages 682–702, https://doi.org/10.1093/gji/ggz423

de Montserrat, A., Aellig, P. S., Schuler, C., Navarrete, I., Räss, L., Fuchs, L., ... & Dominguez, H. (2026). JustRelax. jl: A Julia package for geodynamic modeling with matrix-free solvers. Journal of Open Source Software, 11(118), 9365.

May, D. A., Brown, J., & Le Pourhiet, L. (2015). A scalable, matrix-free multigrid preconditioner for finite element discretizations of heterogeneous Stokes flow. Computer methods in applied mechanics and engineering, 290, 496-523.

Ilangovan, P., Kohl, N., & Mohr, M. (2026). Highly scalable geodynamic simulations with HyTeG. Geoscientific Model Development, 19(4), 1455-1472.

Clevenger, T. C., & Heister, T. (2021). Comparison between algebraic and matrix‐free geometric multigrid for a Stokes problem on adaptive meshes with variable viscosity. Numerical Linear Algebra with Applications, 28(5), e2375.

Reply

Citation: https://doi.org/10.5194/egusphere-2026-2922-RC1

Chase J. Shyu, Sungho Lee, Xuesong Ding, Jaeyoon Keum, Eh Tan, Luc L. Lavier, and Eunseo Choi

Data sets

DynEarthSol (benchmarks-v2.0.0) L. L. Lavier et al. https://doi.org/10.5281/zenodo.20298087

Model code and software

DynEarthSol (v2.0.0) L. L. Lavier et al. https://doi.org/10.5281/zenodo.20293558

Chase J. Shyu, Sungho Lee, Xuesong Ding, Jaeyoon Keum, Eh Tan, Luc L. Lavier, and Eunseo Choi

Viewed

Total article views: 109 (including HTML, PDF, and XML)

HTML	PDF	XML	Total	BibTeX	EndNote
77	27	5	109	4	3

HTML: 77
PDF: 27
XML: 5
Total: 109
BibTeX: 4
EndNote: 3

Views and downloads (calculated since 12 Jun 2026)

Month	HTML	PDF	XML	Total
Jun 2026	47	10	2	59
Jul 2026	30	17	3	50

Cumulative views and downloads (calculated since 12 Jun 2026)

Month	HTML	PDF	XML	Total
Jun 2026	47	10	2	59
Jul 2026	30	17	3	50

Viewed (geographical distribution)

Total article views: 92 (including HTML, PDF, and XML) Thereof 92 with geography defined and 0 with unknown origin.

Country	#	Views	%

Latest update: 24 Jul 2026

Short summary

This study introduces DynEarthSol v2.0, an open-source geoscientific model for simulating tectonic deformation, landscape evolution, and earthquake-cycle behavior. The new version improves computational performance and reproducibility and provides benchmark examples that help other researchers test, reproduce, and extend the model.


Total:	0
HTML:	0
PDF:	0
XML:	0