diff --git a/.github/workflows/draft-pdf.yml b/.github/workflows/draft-pdf.yml
new file mode 100644
index 00000000..15aba992
--- /dev/null
+++ b/.github/workflows/draft-pdf.yml
@@ -0,0 +1,23 @@
+on: [push]
+
+jobs:
+ paper:
+ runs-on: ubuntu-latest
+ name: Paper Draft
+ steps:
+ - name: Checkout
+ uses: actions/checkout@v4
+ - name: Build draft PDF
+ uses: openjournals/openjournals-draft-action@master
+ with:
+ journal: joss
+ # This should be the path to the paper within your repo.
+ paper-path: docs/paper/paper.md
+ - name: Upload
+ uses: actions/upload-artifact@v4
+ with:
+ name: paper
+ # This is the output path where Pandoc will write the compiled
+ # PDF. Note, this should be the same directory as the input
+ # paper.md
+ path: docs/paper/paper.pdf
\ No newline at end of file
diff --git a/.typos.toml b/.typos.toml
index 5e69ca29..f3194677 100644
--- a/.typos.toml
+++ b/.typos.toml
@@ -7,4 +7,4 @@ iy = "iy"
nin = "nin"
[files]
-extend-exclude = ["tutorials/*.pvsm"]
+extend-exclude = ["tutorials/*.pvsm","docs/paper/paper.bib"]
\ No newline at end of file
diff --git a/AUTHORS.md b/AUTHORS.md
index d7277427..14e5df21 100644
--- a/AUTHORS.md
+++ b/AUTHORS.md
@@ -14,8 +14,7 @@ provided substantial additions or modifications. Together, these two groups form
## Contributors
-The following people contributed major additions or modifications to `GeophysicalModelGenerator.jl` and
-are listed in alphabetical order:
+The following people contributed major additions or modifications to `GeophysicalModelGenerator.jl` and are listed in alphabetical order:
* Pascal Aellig
* Albert De Montserrat
diff --git a/Project.toml b/Project.toml
index 4a796a36..7a088ffd 100644
--- a/Project.toml
+++ b/Project.toml
@@ -40,8 +40,8 @@ GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8"
[extensions]
Chmy_utils = "Chmy"
GLMakie_Visualisation = "GLMakie"
-GMT_utils = "GMT"
Gmsh_utils = "GridapGmsh"
+GMT_utils = "GMT"
[compat]
Chmy = "0.1.20"
@@ -51,19 +51,19 @@ FFMPEG = "0.4"
FileIO = "1"
GDAL_jll = "300.900.0 - 301.901.0"
GLMakie = "0.8, 0.9, 0.10"
-GMT = "1.0 - 1.14"
GeoParams = "0.2 - 0.6"
Geodesy = "1"
GeometryBasics = "0.1 - 0.4"
Glob = "1.2 - 1.3"
+GMT = "1"
GridapGmsh = "0.5 - 0.7"
ImageIO = "0.1 - 0.6"
Interpolations = "0.14, 0.15"
JLD2 = "0.4, 0.5"
LightXML = "0.8, 0.9"
MeshIO = "0.1 - 0.4"
-NearestNeighbors = "0.2 - 0.4"
NCDatasets = "0.14"
+NearestNeighbors = "0.2 - 0.4"
Parameters = "0.9 - 0.12"
SpecialFunctions = "1.0, 2"
StaticArrays = "1"
@@ -74,10 +74,9 @@ julia = "1.9"
[extras]
Chmy = "33a72cf0-4690-46d7-b987-06506c2248b9"
GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
-GMT = "5752ebe1-31b9-557e-87aa-f909b540aa54"
KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
[targets]
-test = ["Test", "GMT", "StableRNGs", "GridapGmsh", "Chmy","KernelAbstractions"]
+test = ["Test", "GMT", "StableRNGs", "GridapGmsh", "Chmy", "KernelAbstractions"]
diff --git a/docs/paper/Basic_Tutorial.png b/docs/paper/Basic_Tutorial.png
new file mode 100644
index 00000000..12908952
Binary files /dev/null and b/docs/paper/Basic_Tutorial.png differ
diff --git a/docs/paper/paper.bib b/docs/paper/paper.bib
new file mode 100644
index 00000000..8c89e633
--- /dev/null
+++ b/docs/paper/paper.bib
@@ -0,0 +1,109 @@
+@Article{se-10-1785-2019,
+AUTHOR = {Fraters, M. and Thieulot, C. and van den Berg, A. and Spakman, W.},
+TITLE = {The {G}eodynamic {W}orld {B}uilder: a solution for complex initial conditions in numerical modeling},
+JOURNAL = {Solid Earth},
+VOLUME = {10},
+YEAR = {2019},
+NUMBER = {5},
+PAGES = {1785--1807},
+URL = {https://se.copernicus.org/articles/10/1785/2019/},
+DOI = {10.5194/se-10-1785-2019}
+}
+
+
+ @article{Wessel_Luis_Uieda_Scharroo_Wobbe_Smith_Tian_2019,
+ title={The {G}eneric {M}apping {T}ools {V}ersion 6},
+ volume={20},
+ ISSN={1525-2027, 1525-2027},
+ DOI={10.1029/2019GC008515},
+ abstractNote={The {G}eneric {M}apping {T}ools ({GMT}) software is ubiquitous in the Earth and ocean sciences. As a cross‐platform tool producing high‐quality maps and figures, it is used by tens of thousands of scientists around the world. The basic syntax of GMT scripts has evolved very slowly since the 1990s, despite the fact that GMT is generally perceived to have a steep learning curve with many pitfalls for beginners and experienced users alike. Reducing these pitfalls means changing the interface, which would break compatibility with thousands of existing scripts. With the latest GMT version 6, we solve this conundrum by introducing a new “modern mode” to complement the interface used in previous versions, which GMT 6 now calls “classic mode.” GMT 6 defaults to classic mode and thus is a recommended upgrade for all GMT 5 users. Nonetheless, new users should take advantage of modern mode to make shorter scripts, quickly access commonly used global data sets, and take full advantage of the new tools to draw subplots, place insets, and create animations. Plain Language Summary The Generic Mapping Tools software is widely used in Earth and ocean sciences to process data and make maps and illustrations. This new version simplifies usage, adds quick access to key data sets, and provides a tool for making scientific animations.},
+ number={11},
+ journal={Geochemistry, Geophysics, Geosystems},
+ author={Wessel, P. and Luis, J. F. and Uieda, L. and Scharroo, R. and Wobbe, F. and Smith, W. H. F. and Tian, D.},
+ year={2019}, month=nov, pages={5556–5564}, language={en} }
+
+ @article{DeLaVarga_Schaaf_Wellmann_2019, title={GemPy 1.0: open-source stochastic geological modeling and inversion}, volume={12}, ISSN={1991-9603}, DOI={10.5194/gmd-12-1-2019}, abstractNote={The representation of subsurface structures is an essential aspect of a wide variety of geoscientific investigations and applications, ranging from geofluid reservoir studies, over raw material investigations, to geosequestration, as well as many branches of geoscientific research and applications in geological surveys. A wide range of methods exist to generate geological models. However, the powerful methods are behind a paywall in expensive commercial packages. We present here a full open-source geomodeling method, based on an implicit potential-field interpolation approach. The interpolation algorithm is comparable to implementations in commercial packages and capable of constructing complex full 3-D geological models, including fault networks, fault–surface interactions, unconformities and dome structures. This algorithm is implemented in the programming language Python, making use of a highly efficient underlying library for efficient code generation (Theano) that enables a direct execution on GPUs. The functionality can be separated into the core aspects required to generate 3-D geological models and additional assets for advanced scientific investigations. These assets provide the full power behind our approach, as they enable the link to machine-learning and Bayesian inference frameworks and thus a path to stochastic geological modeling and inversions. In addition, we provide methods to analyze model topology and to compute gravity fields on the basis of the geological models and assigned density values. In summary, we provide a basis for open scientific research using geological models, with the aim to foster reproducible research in the field of geomodeling.}, number={1}, journal={Geoscientific Model Development}, author={De La Varga, Miguel and Schaaf, Alexander and Wellmann, Florian}, year={2019}, month=jan, pages={1–32}, language={en} }
+
+@article{Paffrath_Friederich_Schmid_Handy_2021,
+ title={Imaging structure and geometry of slabs in the greater {A}lpine area – a {P}-wave travel-time tomography using {A}lp{A}rray {S}eismic {N}etwork data},
+ volume={12}, ISSN={1869-9529}, DOI={10.5194/se-12-2671-2021},
+ abstractNote={We perform a teleseismic P-wave travel-time tomography to examine the geometry and structure of subducted lithosphere in the upper mantle beneath the Alpine orogen. The tomography is based on waveforms recorded at over 600 temporary and permanent broadband stations of the dense AlpArray Seismic Network deployed by 24 different European institutions in the greater Alpine region, reaching from the Massif Central to the Pannonian Basin and from the Po Plain to the river Main.}, number={11}, journal={Solid Earth}, author={Paffrath, Marcel and Friederich, Wolfgang and Schmid, Stefan M. and Handy, Mark R. and the AlpArray and AlpArray-Swath D Working Group},
+ year={2021}, month=nov, pages={2671–2702}, language={en} }
+
+ @article{Macherel_Räss_Schmalholz_2024, title={3D Stresses and Velocities Caused by Continental Plateaus: Scaling Analysis and Numerical Calculations With Application to the {T}ibetan Plateau}, volume={25}, ISSN={1525-2027, 1525-2027}, DOI={10.1029/2023GC011356}, abstractNote={Understanding stresses is crucial for geodynamics since they govern rock deformation and metamorphic reactions. However, the magnitudes and distribution of crustal stresses are still uncertain. Here, we use a 3D numerical model in spherical coordinates to investigate stresses and velocities that result from lateral crustal thickness variations around continental plateaus like those observed for the Tibetan plateau. We do not consider any far‐field deformation so that the plateau deforms by horizontal dilatation and vertical thinning. We assume viscous creep, a simplified plateau geometry, and simplified viscosity and density distributions to couple the numerical results with a scaling analysis. Specifically, we study the impact of the viscosity ratio between crust and lithospheric mantle, a rectangular plateau corner, a stress‐dependent power‐law flow law and Earth’s double curvature on the crustal stress field and horizontal velocities. Two orders of magnitude variation in crustal and lithospheric mantle viscosities change the maximum crustal differential stress only by a factor of ≈2. We derive simple analytical estimates for the crustal deviatoric stress and horizontal velocity which agree to first order with 3D numerical calculations. We apply these estimates to calculate the average crustal viscosity in the eastern Tibetan plateau as ≈1022 Pa · s. Furthermore, our results show that a corner strongly affects the stress distribution, particularly the shear stresses, while Earth’s curvature has a minor impact on the stresses. We further discuss potential implications of our results to strike‐slip faulting and fast exhumation around the Tibetan plateau’s syntaxes.}, number={3}, journal={Geochemistry, Geophysics, Geosystems}, author={Macherel, Emilie and Räss, Ludovic and Schmalholz, Stefan M.}, year={2024}, month=mar, pages={e2023GC011356}, language={en} }
+
+ @article{Rappisi_VanderBeek_Faccenda_Morelli_Molinari_2022, title={Slab Geometry and Upper Mantle Flow Patterns in the {C}entral {M}editerranean From {3D} Anisotropic {P}‐Wave Tomography}, volume={127}, ISSN={2169-9313, 2169-9356}, DOI={10.1029/2021JB023488}, abstractNote={We present the first three-dimensional (3D) anisotropic teleseismic P-wave tomography model of the upper mantle covering the entire Central Mediterranean. Compared to isotropic tomography, it is found that including the magnitude, azimuth, and, importantly, dip of seismic anisotropy in our inversions simplifies isotropic heterogeneity by reducing the magnitude of slow anomalies while yielding anisotropy patterns that are consistent with regional tectonics. The isotropic component of our preferred tomography model is dominated by numerous fast anomalies associated with retreating, stagnant, and detached slab segments. In contrast, relatively slower mantle structure is related to slab windows and the opening of back-arc basins. To better understand the complexities in slab geometry and their relationship to surface geological phenomenon, we present a 3D reconstruction of the main Central Mediterranean slabs down to 700 km based on our anisotropic model. P-wave seismic anisotropy is widespread in the Central Mediterranean upper mantle and is strongest at 200–300 km depth. The anisotropy patterns are interpreted as the result of asthenospheric material flowing primarily horizontally around the main slabs in response to pressure exerted by their mid-to-late Cenezoic horizontal motion, while sub-vertical anisotropy possibly reflects asthenospheric entrainment by descending lithosphere. Our results highlight the importance of anisotropic P-wave imaging for better constraining regional upper mantle geodynamics.}, number={5}, journal={Journal of Geophysical Research: Solid Earth}, author={Rappisi, F. and VanderBeek, B. P. and Faccenda, M. and Morelli, A. and Molinari, I.}, year={2022}, month=may, pages={e2021JB023488}, language={en} }
+
+ @phdthesis{Schori_2021, type={PhD}, title={The {D}evelopment of the {J}ura {F}old-and-{T}hrust {B}elt: pre-existing {B}asement {S}tructures and the {F}ormation of {R}amps}, url={https://folia.unifr.ch/unifr/documents/313053},
+ DOI={10.51363/unifr.sth.2022.001},school={University of Fribourg (Switzerland)}, author={Schori, Marc}, year={2021}, month=oct, language={en} }
+
+ @article{McKenzie_1969, title={Speculations on the Consequences and Causes of Plate Motions}, volume={18}, ISSN={0956-540X, 1365-246X}, DOI={10.1111/j.1365-246X.1969.tb00259.x}, abstractNote={Plate theory has successfully related sea floor spreading to the focal mechanisms of earthquakes and the deep structure of island arcs. It is used here to calculate the temperature distribution in the lithosphere thrust beneath island arcs, and to determine the flow and the stress elsewhere in the mantle. Comparison with observations demonstrates that earthquakes are restricted to those regions of the mantle which are colder than a definite temperature. The flow and the stress heating in the mantle can maintain the high heat flow anomaly observed behind island arcs.}, number={1}, journal={Geophysical Journal International}, author={McKenzie, D. P.}, year={1969}, month=sep, pages={1–32}, language={en} }
+
+@article{bezanson2017julia,
+ title={Julia: {A} Fresh Approach to Numerical Computing},
+ author={Bezanson, Jeff and Edelman, Alan and Karpinski, Stefan and
+ Shah, Viral B},
+ journal={SIAM Review},
+ volume={59},
+ number={1},
+ pages={65--98},
+ year={2017},
+ publisher={SIAM},
+ eprint={1411.1607},
+ eprinttype={arxiv},
+ eprintclass={cs.MS},
+ doi={10.1137/141000671}
+}
+
+ @article{desiena2024geophysical,
+ title={Geophysical responses to an environmentally-boosted volcanic unrest},
+ author={De Siena, L and Amoruso, A and Petrosino, S and Crescentini, L},
+ journal={Geophysical Research Letters},
+ volume={51},
+ number={5},
+ pages={e2023GL104895},
+ year={2024},
+ doi={10.1029/2023GL104895},
+ publisher={Wiley Online Library}
+}
+
+@article{gabrielli2023scattering,
+ title={Scattering attenuation images of the control of thrusts and fluid overpressure on the 2016--2017 Central Italy seismic sequence},
+ author={Gabrielli, Simona and Akinci, Aybige and De Siena, Luca and Del Pezzo, Edoardo and Buttinelli, Mauro and Maesano, Francesco Emanuele and Maffucci, Roberta},
+ journal={Geophysical Research Letters},
+ volume={50},
+ number={8},
+ pages={e2023GL103132},
+ year={2023},
+ doi={10.1029/2023GL103132},
+ publisher={Wiley Online Library}
+}
+
+ @inproceedings{Kaus_Popov_Baumann_Pusok_Bauville_Fernandez_Collignon_2016,
+ title={Forward and inverse modelling of lithospheric deformation on geological timescales},
+ volume={48},
+ rights={All rights reserved},
+ ISBN={978-3-95806-109-5},
+ url={https://juser.fz-juelich.de/record/507751/files/nic_2016_kaus.pdf},
+ booktitle={NIC Series},
+ author={Kaus, Boris JP and Popov, Anton A. and Baumann, T. and Püsök, A. and Bauville, Arthur and Fernandez, Naiara and Collignon, Marine},
+ year={2016},
+ pages={299-306}
+}
+@article{napolitano2023imaging,
+ title={Imaging overpressurised fracture networks and geological barriers hindering fluid migrations across a slow-deformation seismic gap},
+ author={Napolitano, Ferdinando and Gabrielli, Simona and De Siena, Luca and Amoroso, Ortensia and Capuano, Paolo},
+ journal={Scientific Reports},
+ volume={13},
+ number={1},
+ pages={19680},
+ year={2023},
+ doi={10.1038/s41598-023-47104-w},
+ publisher={Nature Publishing Group UK London}
+}
+
+ @article{Bauville_Baumann_2019, title={geomIO: An Open‐Source {MATLAB} Toolbox to Create the Initial Configuration of {2‐D/3‐D} Thermo‐Mechanical Simulations From {2‐D} Vector Drawings}, volume={20}, ISSN={1525-2027, 1525-2027}, DOI={10.1029/2018GC008057}, number={3}, journal={Geochemistry, Geophysics, Geosystems}, author={Bauville, A. and Baumann, T. S.}, year={2019}, month=mar, pages={1665–1675}, language={en} }
+
+ @article{Shephard_Matthews_Hosseini_Domeier_2017, title={On the consistency of seismically imaged lower mantle slabs}, volume={7}, ISSN={2045-2322}, DOI={10.1038/s41598-017-11039-w}, abstractNote={Abstract
+ The geoscience community is increasingly utilizing seismic tomography to interpret mantle heterogeneity and its links to past tectonic and geodynamic processes. To assess the robustness and distribution of positive seismic anomalies, inferred as subducted slabs, we create a set of vote maps for the lower mantle with 14 global P-wave or S-wave tomography models. Based on a depth-dependent threshold metric, an average of 20% of any given tomography model depth is identified as a potential slab. However, upon combining the 14 models, the most consistent positive wavespeed features are identified by an increasing vote count. An overall peak in the most robust anomalies is found between 1000–1400 km depth, followed by a decline to a minimum around 2000 km. While this trend could reflect reduced tomographic resolution in the middle mantle, we show that it may alternatively relate to real changes in the time-dependent subduction flux and/or a mid-lower mantle viscosity increase. An apparent secondary peak in agreement below 2500 km depth may reflect the degree-two lower mantle slow seismic structures. Vote maps illustrate the potential shortcomings of using a limited number or type of tomography models and slab threshold criteria.}, number={1}, journal={Scientific Reports}, author={Shephard, G. E. and Matthews, K. J. and Hosseini, K. and Domeier, M.}, year={2017}, month=sep, pages={10976}, language={en} }
+
diff --git a/docs/paper/paper.md b/docs/paper/paper.md
index b250d2b3..ec15e136 100644
--- a/docs/paper/paper.md
+++ b/docs/paper/paper.md
@@ -1,11 +1,12 @@
---
-title: 'GeophysicalModelGenerator.jl: A Julia package to visualize geoscientific data and generate model '
+title: 'GeophysicalModelGenerator.jl: A Julia package to visualise geoscientific data and create numerical model setups'
tags:
- - Julia
+ - julia
- geosciences
- geodynamics
- tectonics
- geophysics
+ - computational geosciences
authors:
- name: Boris J.P. Kaus
orcid: 0000-0002-0247-8660
@@ -13,107 +14,200 @@ authors:
- name: Marcel Thielmann
orcid: 0000-0003-1185-3730
affiliation: 2
- - name: Arne Spang
- orcid: 0000-0002-6107-0403
- affiliation: 2
- - name: Luca de Siena
- orcid: 0000-0002-3615-5923
- affiliation: 3
- name: Pascal Aellig
orcid: 0009-0008-9039-5646
affiliation: 1
+ - name: Albert de Montserrat
+ orcid: 0000-0003-1694-3735
+ affiliation: 3
+ - name: Luca de Siena
+ orcid: 0000-0002-3615-5923
+ affiliation: 4,1
- name: Jacob Frasukiewicz
- #orcid: 0000-0003-1185-3730
+ orcid: 0009-0002-5049-4259
+ affiliation: 1
+ - name: Lukas Fuchs
+ orcid: 0000-0002-9165-6384
+ affiliation: 5
+ - name: Andrea Piccolo
+ orcid: 0000-0003-3074-6041
+ affiliation: 2
+ - name: Hendrik Ranocha
+ orcid: 0000-0002-3456-2277
affiliation: 1
+ - name: Nicolas Riel
+ orcid: 0000-0002-5037-5519
+ affiliation: 1
+ - name: Christian Schuler
+ orcid: 0009-0004-9873-9774
+ affiliation: 1
+ - name: Arne Spang
+ orcid: 0000-0002-6107-0403
+ affiliation: 2
+ - name: Tatjana Weiler
+ #orcid: 0000-0002-6107-0403
+ affiliation: 2
+
affiliations:
- - name: Johannes Gutenberg University Mainz, Germany
+ - name: Johannes Gutenberg-University Mainz, Germany
index: 1
- name: University of Bayreuth, Germany
index: 2
- - name: Bologna University, Italy
+ - name: ETH Zürich, Switzerland
index: 3
+ - name: Alma Mater Studiorum Bologna University, Italy
+ index: 4
+ - name: Goethe University Frankfurt, Germany
+ index: 5
-date: 27 February 2024
+date: 14 March 2024
bibliography: paper.bib
---
# Summary
-The forces on stars, galaxies, and dark matter under external gravitational
-fields lead to the dynamical evolution of structures in the universe. The orbits
-of these bodies are therefore key to understanding the formation, history, and
-future state of galaxies. The field of "galactic dynamics," which aims to model
-the gravitating components of galaxies to study their structure and evolution,
-is now well-established, commonly taught, and frequently used in astronomy.
-Aside from toy problems and demonstrations, the majority of problems require
-efficient numerical tools, many of which require the same base code (e.g., for
-performing numerical orbit integration).
+Geoscientific data exists in a wide variety of formats. Yet, to make a consistent interpretation of a certain region, it is often helpful to jointly visualise all this data using the same coordinates and compare, for example, seismic tomography, surface geology, Moho (crust-mantle transition) depth, Earthquake locations, and GPS surface velocities. If one wishes to create mechanical or thermo-mechanical numerical models of the region, creating an input setup that honors these constraints is crucial. And since most numerical codes work in Cartesian boxes, it is helpful to have tools to project the data from geographic to Cartesian coordinates.
-# Statement of need
+A significant challenge in doing this is that there is no standard format for geoscientific data. Seismic tomography, for example, may come in the form of ASCII data with `lon/lat/depth` axes or as NetCDF files, with the ordering of the data typically differing from one dataset to the other. In ideal cases, geological surfaces may be provided as GeoTIFF images. In many cases, however, the underlying data discussed in publications are not available in digital format and are only shown as figures in the paper. It is nevertheless still helpful to visualise these in 3D in the correct coordinates, along with more recent, digitally available datasets.
-`Gala` is an Astropy-affiliated Python package for galactic dynamics. Python
-enables wrapping low-level languages (e.g., C) for speed without losing
-flexibility or ease-of-use in the user-interface. The API for `Gala` was
-designed to provide a class-based and user-friendly interface to fast (C or
-Cython-optimized) implementations of common operations such as gravitational
-potential and force evaluation, orbit integration, dynamical transformations,
-and chaos indicators for nonlinear dynamics. `Gala` also relies heavily on and
-interfaces well with the implementations of physical units and astronomical
-coordinate systems in the `Astropy` package [@astropy] (`astropy.units` and
-`astropy.coordinates`).
+The aim of the `GeophysicalModelGenerator.jl` package is therefore twofold:
-`Gala` was designed to be used by both astronomical researchers and by
-students in courses on gravitational dynamics or astronomy. It has already been
-used in a number of scientific publications [@Pearson:2017] and has also been
-used in graduate courses on Galactic dynamics to, e.g., provide interactive
-visualizations of textbook material [@Binney:2008]. The combination of speed,
-design, and support for Astropy functionality in `Gala` will enable exciting
-scientific explorations of forthcoming data releases from the *Gaia* mission
-[@gaia] by students and experts alike.
+1) Simplify collecting and visualising a wide variety of geoscientific data that is provided as point (e.g., earthquake locations), surface (e.g., topography) or volumetric data (e.g., seismic tomography).
+2) Create input setups for 2D or 3D numerical models.
-# Mathematics
+# Statement of need
-Single dollars ($) are required for inline mathematics e.g. $f(x) = e^{\pi/x}$
+`GeophysicalModelGenerator.jl` is a Julia [@bezanson2017julia] package that helps collect and visualise a wide variety of geophysical and geoscientific data coherently. It also simplifies the process of generating 2D or 3D models that can, for example, be used as input models in geodynamic simulations. It provides functions that transfer data from one format to the other, or project them from geographic `Longitude/Latitude/Depth` or `UTM` coordinates to Cartesian coordinates (`kilometers`). It allows performing tasks such as creating cross-sections through volumetric data, importing screenshots from published papers, downloading digital elevation data and saving the resulting data in `VTK` format, which can, for instance, be visualised with open source tools such as [Paraview](www.paraview.org).
+
+Most geoscientists tend to have their own Python/MATLAB/Bash visualisation, and many perform part of this job on a daily basis. Yet, having all functionality in one place in an easy-to-use package makes this more extendable and will likely facilitate sharing data along with their interpretations.
+
+# Related software packages
+Perhaps the most widely used package in geophysics to create figures or maps is the Generic Mapping Tools package ([`GMT`](https://www.generic-mapping-tools.org)), which also provides a Julia interface [GMT.jl](https://github.com/GenericMappingTools/GMT.jl) [@Wessel_Luis_Uieda_Scharroo_Wobbe_Smith_Tian_2019]. It mostly focuses on generating maps and postscript/pdf images and is therefore not ideally suited for interactive 3D data visualisation, or to generate input models for numerical codes.
+
+The [`Geodynamic World Builder`](https://github.com/GeodynamicWorldBuilder/WorldBuilder) is a C++ library to create model setups [@se-10-1785-2019]. The focus is on generating input models for geodynamic simulations, such as subduction zones and related thermal structures. It has C and Fortran wrappers and can thus be embedded in geodynamic codes. Users of the `Geodynamic World Builder` have to generate JSON files to define the model geometry, which is less interactive than by using the Julia `REPL`. There is no straightforward way to integrate existing geophysical/geological data in the workflow and compare model results with them.
+
+[geomIO](https://bitbucket.org/geomio/geomio/src) is a MATLAB-based toolbox that allows creating geodynamic input setups by drawing several cross-sections using vector software such as Inkscape, which is put together into 3D volumes [@Bauville_Baumann_2019]. While it does allow the creation of sophisticated setups, data can only be taken into account by adding them as screenshots to Inkscape while drawing cross-sections. Its reliance on commercial software may be problematic for some users.
+
+[GemPy](https://www.gempy.org) is a Python-based, open-source geomodeling library that can construct 3D geological models of folded structures, fault networks and unconformities while taking uncertainties into account [@DeLaVarga_Schaaf_Wellmann_2019]. It focuses on creating geometric models with uncertainties rather than on integrating a wide variety of geoscientific datasets.
+
+There are also a number of commercial software solutions:
+
+- [Petrel subsurface software](https://www.software.slb.com/products/petrel) (by Schlumberger), which is mostly used by the hydrocarbon industry and is particularly powerful in integrating seismic reflection and well-data.
+
+- [GOCAD Mining Suite](https://www.mirageoscience.com/mining-industry-software/gocad-mining-suite/) (by MiraGeoscience) helps to generate geometric models of the subsurface in the vicinity of mines based on sparse geological measurements and drill hole data.
+
+- [GeoModeller](https://www.intrepid-geophysics.com/products/geomodeller/) (by Intrepid Geophysics) creates surface-near geometric geological models by implicit modelling of surface measurements while taking geophysical constraints into account.
+
+In all cases, the commercial license fees are far beyond what most researchers can afford, even if reduced license fees are often available for academia. The closed-source nature of the software packages makes them non-extendable by the community.
+
+The `GeophysicalModelGenerator.jl` package is already used to generate input models for the geodynamic codes [LaMEM](https://github.com/UniMainzGeo/LaMEM) [@Kaus_Popov_Baumann_Pusok_Bauville_Fernandez_Collignon_2016], [JustRelax.jl](https://github.com/PTsolvers/JustRelax.jl), and [MagmaThermokinematics.jl](https://github.com/boriskaus/MagmaThermoKinematics.jl). It is also used in a number of short courses and lectures at the Universities of Mainz, Heidelberg, Jena and Bologna, in hands-on tutorials during various workshops and in a number of recent publications to highlight 3D [@gabrielli2023scattering; @napolitano2023imaging] and 4D [@desiena2024geophysical] geophysical and geological data.
-Double dollars make self-standing equations:
+# Basic usage
+
+The core of the package consists of the `GeoData`, `UTMData`, `ParaviewData`, and `CartData` structures which hold the 3D data along with coordinates (and potentially metadata) information. `GeophysicalModelGenerator.jl` can be installed using the built-in Julia package manager:
+```julia
+julia> ]
+(@v1.10) pkg> add GeophysicalModelGenerator
+```
+which comes with a test-suite:
+```julia
+(@v1.10) pkg> test GeophysicalModelGenerator
+```
+and can be loaded with:
+```julia
+julia> using GeophysicalModelGenerator
+```
-$$\Theta(x) = \left\{\begin{array}{l}
-0\textrm{ if } x < 0\cr
-1\textrm{ else}
-\end{array}\right.$$
+As a first example, we will download a 3D seismic tomography dataset for the Alpine region (from [@Paffrath_Friederich_Schmid_Handy_2021]):
+```julia
+julia> Tomo_Alps_full = load_GMG(
+ "https://zenodo.org/records/10738510/files/Paffrath_2021_SE_Pwave.jld2?download=1");
+```
+We can download the topography of the Alpine region with:
+```julia
+julia> Topo_Alps = load_GMG(
+ "https://zenodo.org/records/10738510/files/AlpsTopo.jld2?download=1");
+```
-You can also use plain \LaTeX for equations
-\begin{equation}\label{eq:fourier}
-\hat f(\omega) = \int_{-\infty}^{\infty} f(x) e^{i\omega x} dx
-\end{equation}
-and refer to \autoref{eq:fourier} from text.
+The seismic data covers a much wider region than the Alps itself, but in much of that region there is poor data coverage. We can, therefore, extract a part of the data that has coverage:
+```julia
+julia> Tomo_Alps = extract_subvolume(Tomo_Alps_full, Lon_level=(4,20),
+ Lat_level=(36,50), Depth_level=(-600,-10));
+```
-# Citations
+At this stage, we can save the data to `VTK` format:
+```julia
+julia> write_paraview(Tomo_Alps,"Tomo_Alps");
+julia> write_paraview(Topo_Alps,"Topo_Alps");
+```
+and open it with Paraview (see \autoref{fig:basic}a).
+We can create vertical and horizontal cross-sections through the data with:
+```julia
+julia> Cross_200km = cross_section(Tomo_Alps, Depth_level=-200, Interpolate=true);
+julia> Cross_vert = cross_section(Tomo_Alps, Start=(5,47), End=(15,44));
+julia> write_paraview(Cross_vert, "Cross_vert");
+julia> write_paraview(Cross_200km,"Cross_200km");
+```
+and visualise them along with the volumetric data (\autoref{fig:basic}a).
-Citations to entries in paper.bib should be in
-[rMarkdown](http://rmarkdown.rstudio.com/authoring_bibliographies_and_citations.html)
-format.
+![Example of visualising 3D seismic data of the Alps, using a) geographic coordinates (`GeoData`) or b) Cartesian coordinates (`CartData`) projected from geographic coordinates. The topography and several slices through the 3D seismic tomography P-wave model of [@Paffrath_Friederich_Schmid_Handy_2021] are shown. \label{fig:basic} ](Basic_Tutorial.png){ width=100% }
-If you want to cite a software repository URL (e.g. something on GitHub without a preferred
-citation) then you can do it with the example BibTeX entry below for @fidgit.
+One complication with geographic data is that Paraview does not have native support for geographic coordinates; accordingly, it is not always straightforward to use the built-in tools, for example, to create slices through the data. In addition, many numerical models work in (orthogonal) Cartesian rather than in spherical coordinates, which appears to be a good first-order approximation for many geodynamic applications [@Macherel_Räss_Schmalholz_2024].
-For a quick reference, the following citation commands can be used:
-- `@author:2001` -> "Author et al. (2001)"
-- `[@author:2001]` -> "(Author et al., 2001)"
-- `[@author1:2001; @author2:2001]` -> "(Author1 et al., 2001; Author2 et al., 2002)"
+`GeophysicalModelGenerator.jl` includes tools to transfer the data from geographic to Cartesian coordinates, which requires defining a projection point along which the projection is performed:
+```julia
+julia> proj = ProjectionPoint(Lon=12.0,Lat =43)
+ProjectionPoint(43.0, 12.0, 255466.98055255096, 4.765182932801006e6, 33, true)
+```
+We can now project the topography with:
+```julia
+julia> Topo_cart = convert2CartData(Topo_Alps, proj);
+```
+which returns a `CartData` (Cartesian data) structure. The disadvantage of doing this projection is that the resulting Cartesian grid is no longer strictly orthogonal, which is a problem for some Cartesian numerical models (e.g., those using a finite difference discretisation).
+We can project the data on an orthogonal grid as well, by first creating appropriately sized orthogonal grids for the tomography and topography:
+```julia
+julia> Tomo_rect = CartData(xyz_grid(-550.0:10:600, -500.0:10:700, -600.0:5:-17));
+julia> Topo_rect = CartData(xyz_grid(-550.0:1:600, -500.0:1:700, 0));
+```
+Next, we can project the data to the orthogonal grids with:
+```julia
+julia> Topo_rect = project_CartData(Topo_rect, Topo_Alps, proj);
+julia> Tomo_rect = project_CartData(Tomo_rect, Tomo_Alps, proj);
+julia> write_paraview(Tomo_rect,"Tomo_rect");
+julia> write_paraview(Topo_rect,"Topo_rect");
+```
+We can now use the built-in tools of Paraview to visualise the data (see \autoref{fig:basic} b), and use this as inspiration to create an initial numerical model setup. It is also possible to interpolate other seismic tomography datasets to the same grid and subsequently compute a _votemap_ to count how many tomographic models contain a specific seismic anomaly [see, e.g., @Shephard_Matthews_Hosseini_Domeier_2017].
-# Figures
+# Examples of usage
+`GeophysicalModelGenerator.jl` comes with build-in (CI/CD) tests and [tutorials](https://juliageodynamics.github.io/GeophysicalModelGenerator.jl/dev) that explain the most important use cases, from importing data to generating input model setups for numerical simulations. In the following, we present a number of examples that illustrate various aspects of the package. Many additional tutorials are available in the online documentation.
-Figures can be included like this:
-
-and referenced from text using \autoref{fig:example}.
+### Visualise data from the Alps
+The European Alps are among the best-studied mountain belts on the planet and have therefore been the focus of numerous geological and geophysical studies. Different seismic tomography models have been published (using different parameterisations and datasets), but those do not necessarily agree with each other.
+
+In [`Tutorial_AlpineData.jl`](https://juliageodynamics.github.io/GeophysicalModelGenerator.jl/dev/man/Tutorial_AlpineData/), users learn how to load the topography of the region, import Moho data, load and visualise GPS vectors, import and plot earthquake locations, along with cross-sections through the model (\autoref{fig:alps}).
-Figure sizes can be customized by adding an optional second parameter:
-{ width=20% }
+![Example of combined data of the Alps, which shows the GPS surface velocity (arrows), topography, earthquake locations (colored dots) and cross-sections through a recent anisotropic P-wave tomography model by [@Rappisi_VanderBeek_Faccenda_Morelli_Molinari_2022]. \label{fig:alps}](../src/assets/img/GMG_AlpineData.png){ width=90% }
+
+### La Palma volcanic eruption
+The 2019 Cumbre Viejo eruption in La Palma, Canary Islands, was accompanied by seismic activity. In [`Tutorial_LaPalma.jl`](https://juliageodynamics.github.io/GeophysicalModelGenerator.jl/dev/man/Tutorial_LaPalma/), users learn to generate a Cartesian block model of the island, import seismicity and use that to generate a 3D volumetric seismic activity map (\autoref{fig:lapalma}).
+
+{ width=100% }
+
+
+### Jura mountains
+The Jura mountains are a small-scale fold and thrust belt located in Switzerland and France. Thanks to seismic cross-sections and boreholes, a lot of information is available about its structure at depth. This information was used to generate extensive 3D models of the subsurface, including thickness maps of various geological units, generate a new geological map of the region, and create balanced reconstructions [@Schori_2021].
+
+In [`Tutorial_Jura.jl`](https://juliageodynamics.github.io/GeophysicalModelGenerator.jl/dev/man/Tutorial_Jura/), users learn how to drape the geological map over the topography, import surfaces from GeoTIFF images (such as basement topography), and include screenshots from geological cross-sections. The data is rotated and transferred to Cartesian coordinates such that we obtain a 3D block model that is perpendicular to the strike of the mountain range (\autoref{fig:jura}).
+
+![Example of creating a 3D Cartesian block model that runs perpendicular to the Jura mountains, combining surface geology, with screenshots from interpreted cross-sections (in the center right), and digital data of the the basement topography [using data of @Schori_2021]. \label{fig:jura}](../src/assets/img/Jura_2.png){ width=100% }
+
+
+### Slab model setup
+In `Tutorial_NumericalModel_3D.jl`, users learn how to generate a 3D geodynamic model setup with subducting slabs, a mid-oceanic ridge and an overriding cratonic lithosphere. The thermal structure of the oceanic slab increases away from the ridge until the trench, following a halfspace cooling analytical solution. In contrast, the thermal structure of the subducted part of the slab is based on an analytical solution that takes heating from the surrounding hot mantle into account [after @McKenzie_1969]. Between the mantle and the trench, the slab uses a mixture of these two thermal models. A weak zone is added above the slab (to facilitate subduction in numerical models). A sedimentary wedge is situated at the continental margin, and a grid-like pattern is put on top of the oceanic slab to simplify tracking deformation throughout the simulation (\autoref{fig:slab3d}).
+
+{ width=100% }
# Acknowledgements
-
-We acknowledge contributions from Brigitta Sipocz, Syrtis Major, and Semyeong
-Oh, and support from Kathryn Johnston during the genesis of this project.
+We acknowledge funding by 1) the European Research Council as part of Consolidator Grant 771143 (MAGMA), 2) the German Research Foundation through DFG grants KA3367/10-1, TH2076/7-1 (both part of the SPP 2017 4D-MB project) and Emmy Noether grant TH2076/8-1, 3) the German Ministry of Science and Education (BMBF) through project DEGREE, and 4) by the CHEESE-2p Center of Excellence, co-funded by both EuroHPC-JU and the BMBF.
# References
diff --git a/docs/src/assets/img/Basic_Tutorial.png b/docs/src/assets/img/Basic_Tutorial.png
new file mode 100644
index 00000000..b1129eba
Binary files /dev/null and b/docs/src/assets/img/Basic_Tutorial.png differ
diff --git a/docs/src/assets/img/GMG_AlpineData.png b/docs/src/assets/img/GMG_AlpineData.png
index ee0e4678..1675ff1a 100644
Binary files a/docs/src/assets/img/GMG_AlpineData.png and b/docs/src/assets/img/GMG_AlpineData.png differ
diff --git a/docs/src/assets/img/Jura_2.png b/docs/src/assets/img/Jura_2.png
index 93b5cbe0..8c4bca7f 100644
Binary files a/docs/src/assets/img/Jura_2.png and b/docs/src/assets/img/Jura_2.png differ
diff --git a/docs/src/assets/img/Mechanical3D_Tutorial_2.png b/docs/src/assets/img/Mechanical3D_Tutorial_2.png
index 91395102..b16248ce 100644
Binary files a/docs/src/assets/img/Mechanical3D_Tutorial_2.png and b/docs/src/assets/img/Mechanical3D_Tutorial_2.png differ
diff --git a/docs/src/assets/img/Tutorial_LaPalma.png b/docs/src/assets/img/Tutorial_LaPalma.png
index c3a0bfb7..9cf83d37 100644
Binary files a/docs/src/assets/img/Tutorial_LaPalma.png and b/docs/src/assets/img/Tutorial_LaPalma.png differ
diff --git a/ext/GMT_utils.jl b/ext/GMT_utils.jl
index a502d82b..77878ec8 100644
--- a/ext/GMT_utils.jl
+++ b/ext/GMT_utils.jl
@@ -68,8 +68,8 @@ julia> write_paraview(Topo,"Topo_Alps")
function import_topo(limits; file::String="@earth_relief_01m", maxattempts=5)
# Correct if negative values are given (longitude coordinates that are west)
- ind = limits[1:2] .< 0
-
+ ind = findall(limits[1:2] .< 0);
+
if (limits[1] < 0) && (limits[2] < 0)
limits[ind] .= 360 .+ limits[ind]
limits[1:2] = sort(limits[1:2])
@@ -88,9 +88,7 @@ function import_topo(limits; file::String="@earth_relief_01m", maxattempts=5)
end
attempt += 1
end
- if isempty(G)
- error("Could not download GMT topography data")
- end
+ (@isdefined G) || error("Could not download GMT topography data")
# Transfer to GeoData
nx, ny = size(G.z,2), size(G.z,1)
diff --git a/test/runtests.jl b/test/runtests.jl
index 29b9a822..bcf3ef42 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -87,6 +87,9 @@ using Test
end
end
+# Include tutorials
+include("test_tutorials.jl")
+
# Cleanup
foreach(rm, filter(endswith(".vts"), readdir()))
foreach(rm, filter(endswith(".vtu"), readdir()))
diff --git a/test/test_tutorials.jl b/test/test_tutorials.jl
new file mode 100644
index 00000000..99c16ba4
--- /dev/null
+++ b/test/test_tutorials.jl
@@ -0,0 +1,26 @@
+#this runs the tutorials that are part of the JOSS paper to ensure that they keep working moving forward
+
+@testset "Basic tutorial" begin
+ include("../tutorials/Tutorial_Basic.jl")
+end
+
+#@testset "Jura tutorial" begin
+# include("../tutorials/Tutorial_Jura.jl")
+#end
+
+@testset "LaPalma tutorial" begin
+ include("../tutorials/Tutorial_LaPalma.jl")
+end
+
+# Deactivating this one as it plots in the tutorial
+#@testset "AlpineData tutorial" begin
+# include("../tutorials/Tutorial_AlpineData.jl")
+#end
+
+@testset "2D Numerical Model tutorial" begin
+ include("../tutorials/Tutorial_NumericalModel_2D.jl")
+end
+#
+#@testset "3D Numerical Model tutorial" begin
+# include("../tutorials/Tutorial_NumericalModel_3D.jl")
+#end
\ No newline at end of file
diff --git a/tutorials/Jura_2.pvsm b/tutorials/Jura_2.pvsm
index db1d3907..99a175ab 100644
--- a/tutorials/Jura_2.pvsm
+++ b/tutorials/Jura_2.pvsm
@@ -60,7 +60,7 @@
-
+
@@ -104,48 +104,48 @@
-
-
+
+
-
+
-
+
-
+
-
+
-