Numerical model setups

In order to generate numerical simulations from setups created with GeophysicalModelGenerator.jl, we provide a few routines that directly create setups.

The routines provided here have the following functionality:

Add lithospheric boxes to a setup, that may have a layered structure and various thermal structures
Add various geometries (spheres, cuylinders, ellipsoids)
Add lithospheric structure
Add various 1D thermal structures (and possibilities to combine them)

GeophysicalModelGenerator.add_box! — Function

add_box!(Phase, Temp, Grid::AbstractGeneralGrid; xlim::Tuple = (20,100), [ylim::Tuple = (1,10)], zlim::Tuple = (10,80),
        Origin=nothing, StrikeAngle=0, DipAngle=0,
        phase = ConstantPhase(1),
        T=nothing,
        cell=false )

Adds a box with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models

Parameters

Phase - Phase array (consistent with Grid)
Temp - Temperature array (consistent with Grid)
Grid - grid structure (can be any of the grid types in GMG)
xlim - left/right coordinates of box
ylim - front/back coordinates of box [optional; if not specified we use the whole box]
zlim - bottom/top coordinates of box
Origin - the origin, used to rotate the box around. Default is the left-front-top corner
StrikeAngle - strike angle of slab
DipAngle - dip angle of slab
phase - specifies the phase of the box. See ConstantPhase(),LithosphericPhases()
T - specifies the temperature of the box. See ConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp(),LithosphericTemp()
cell - if true, Phase and Temp are defined on centers

Examples

Example 1) Box with constant phase and temperature & a dip angle of 10 degrees:

julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> add_box!(Phases,Temp,Grid, xlim=(0,500), zlim=(-50,0), phase=ConstantPhase(3), DipAngle=10, T=ConstantTemp(1000))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

Example 2) Box with halfspace cooling profile

julia> Grid = CartData(xyz_grid(-1000:10:1000,0,-660:10:0))
julia> Phases = zeros(Int32,   size(Grid));
julia> Temp   = zeros(Float64, size(Grid));
julia> add_box!(Phases,Temp,Grid, xlim=(0,500), zlim=(-50,0), phase=ConstantPhase(3), DipAngle=10, T=HalfspaceCoolingTemp(Age=30))
julia> Grid = addfield(Grid, (;Phases,Temp));       # Add to Cartesian model
julia> write_paraview(Grid,"LaMEM_ModelSetup")  # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

source

GeophysicalModelGenerator.add_layer! — Function

add_layer!(Phase, Temp, Grid::AbstractGeneralGrid; xlim::Tuple = (1,100), [ylim::Tuple = (0,20)], zlim::Tuple = (0,-100),
        phase = ConstantPhase(1),
        T=nothing, cell=false )

Adds a layer with phase & temperature structure to a 3D model setup. The most common use would be to add a lithospheric layer to a model setup. This simplifies creating model geometries in geodynamic models

Parameters

Phase - Phase array (consistent with Grid)
Temp - Temperature array (consistent with Grid)
Grid - grid structure (usually obtained with readLaMEMinputfile, but can also be other grid types)
xlim - left/right coordinates of box
ylim - front/back coordinates of box
zlim - bottom/top coordinates of box
phase - specifies the phase of the box. See ConstantPhase(),LithosphericPhases()
T - specifies the temperature of the box. See ConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp()

Examples

Example 1) Layer with constant phase and temperature

julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> add_layer!(Phases,Temp,Grid, zlim=(-50,0), phase=ConstantPhase(3), T=ConstantTemp(1000))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

Example 2) Box with halfspace cooling profile

julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> add_layer!(Phases,Temp,Grid, zlim=(-50,0), phase=ConstantPhase(3), T=HalfspaceCoolingTemp())
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

source

GeophysicalModelGenerator.add_sphere! — Function

add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid; cen::Tuple = (0,0,-1), radius::Number,
        phase = ConstantPhase(1).
        T=nothing, cell=false )

Adds a sphere with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models

Parameters

Phase - Phase array (consistent with Grid)
Temp - Temperature array (consistent with Grid)
Grid - LaMEM grid structure (usually obtained with readLaMEMinputfile)
cen - center coordinates of sphere
radius - radius of sphere
phase - specifies the phase of the box. See ConstantPhase(),LithosphericPhases()
T - specifies the temperature of the box. See ConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp()
cell - if true, Phase and Temp are defined on cell centers

Example

Sphere with constant phase and temperature:

julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> add_sphere!(Phases,Temp,Grid, cen=(0,0,-1), radius=0.5, phase=ConstantPhase(2), T=ConstantTemp(800))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

source

GeophysicalModelGenerator.add_ellipsoid! — Function

add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; cen::Tuple = (-1,-1,-1), axes::Tuple = (0.2,0.1,0.5),
        Origin=nothing, StrikeAngle=0, DipAngle=0,
        phase = ConstantPhase(1).
        T=nothing, cell=false )

Adds an Ellipsoid with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models

Parameters

Phase - Phase array (consistent with Grid)
Temp - Temperature array (consistent with Grid)
Grid - LaMEM grid structure (usually obtained with readLaMEMinputfile)
cen - center coordinates of sphere
axes - semi-axes of ellipsoid in X,Y,Z
Origin - the origin, used to rotate the box around. Default is the left-front-top corner
StrikeAngle - strike angle of slab
DipAngle - dip angle of slab
phase - specifies the phase of the box. See ConstantPhase(),LithosphericPhases()
T - specifies the temperature of the box. See ConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp()
cell - if true, Phase and Temp are defined on cell centers

Example

Ellipsoid with constant phase and temperature, rotated 90 degrees and tilted by 45 degrees:

julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> add_ellipsoid!(Phases,Temp,Grid, cen=(-1,-1,-1), axes=(0.2,0.1,0.5), StrikeAngle=90, DipAngle=45, phase=ConstantPhase(3), T=ConstantTemp(600))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

source

GeophysicalModelGenerator.add_cylinder! — Function

add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; base::Tuple = (-1,-1,-1.5), cap::Tuple = (-1,-1,-0.5), radius::Number,
        phase = ConstantPhase(1),
        T=nothing, cell=false )

Adds a cylinder with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models

Parameters

Phase - Phase array (consistent with Grid)
Temp - Temperature array (consistent with Grid)
Grid - Grid structure (usually obtained with read_LaMEM_inputfile)
base - center coordinate of bottom of cylinder
cap - center coordinate of top of cylinder
radius - radius of the cylinder
phase - specifies the phase of the box. See ConstantPhase(),LithosphericPhases()
T - specifies the temperature of the box. See ConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp()
cell - if true, Phase and Temp are defined on cell centers

Example

Cylinder with constant phase and temperature:

julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> add_cylinder!(Phases,Temp,Grid, base=(-1,-1,-1.5), cap=(1,1,-0.5), radius=0.25, phase=ConstantPhase(4), T=ConstantTemp(400))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

source

GeophysicalModelGenerator.add_stripes! — Function

add_stripes!(Phase, Grid::AbstractGeneralGrid;
    stripAxes       = (1,1,0),
    stripeWidth     =  0.2,
    stripeSpacing   =  1,
    Origin          =  nothing,
    StrikeAngle     =  0,
    DipAngle        =  10,
    phase           =  ConstantPhase(3),
    stripePhase     =  ConstantPhase(4),
    cell            = false)

Adds stripes to a pre-defined phase (e.g. added using add_box!)

Parameters

Phase - Phase array (consistent with Grid)
Grid - grid structure (usually obtained with readLaMEMinputfile, but can also be other grid types)
stripAxes - sets the axis for which we want the stripes. Default is (1,1,0) i.e. X, Y and not Z
stripeWidth - width of the stripe
stripeSpacing - space between two stripes
Origin - the origin, used to rotate the box around. Default is the left-front-top corner
StrikeAngle - strike angle
DipAngle - dip angle
phase - specifies the phase we want to apply stripes to
stripePhase - specifies the stripe phase
cell - if true, Phase and Temp are defined on centers

Example

Example: Box with striped phase and constant temperature & a dip angle of 10 degrees:

julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> add_box!(Phases,Temp,Grid, xlim=(0,500), zlim=(-50,0), phase=ConstantPhase(3), DipAngle=10, T=ConstantTemp(1000))
julia> add_stripes!(Phases, Grid, stripAxes=(1,1,1), stripeWidth=0.2, stripeSpacing=1, Origin=nothing, StrikeAngle=0, DipAngle=10, phase=ConstantPhase(3), stripePhase=ConstantPhase(4))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

source

GeophysicalModelGenerator.add_slab! — Function

add_slab!(Phase, Temp, Grid::AbstractGeneralGrid,  trench::Trench; phase = ConstantPhase(1), T = nothing, cell=false)

Adds a curved slab with phase & temperature structure to a 3D model setup.

Parameters

Phase - Phase array (consistent with Grid)
Temp - Temperature array (consistent with Grid)
Grid - grid structure (can be any of the grid types in GMG)
trench - Trench structure
phase - specifies the phase of the box. See ConstantPhase(),LithosphericPhases()
T - specifies the temperature of the box. See ConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp(),LithosphericTemp()
cell - if true, Phase and Temp are defined on cells

Examples

Example 1) Slab

julia> x     = LinRange(0.0,1200.0,128);
julia> y     = LinRange(0.0,1200.0,128);
julia> z     = LinRange(-660,50,128);
julia> Cart  = CartData(xyz_grid(x, y, z));
julia> Phase = ones(Int64,size(Cart));
julia> Temp  = fill(1350.0,size(Cart));
# Define the trench:
julia> trench= Trench(Start = (400.0,400.0), End = (800.0,800.0), θ_max = 45.0, direction = 1.0, n_seg = 50, Length = 600.0, Thickness = 80.0, Lb = 500.0, d_decoupling = 100.0, type_bending =:Ribe)
julia> phase = LithosphericPhases(Layers=[5 7 88], Phases = [2 3 4], Tlab=nothing)
julia> TsHC  = HalfspaceCoolingTemp(Tsurface=20.0, Tmantle=1350, Age=30, Adiabat=0.4)
julia> add_slab!(Phase, Temp, Cart, trench, phase = phase, T = TsHC)

source

GeophysicalModelGenerator.make_volc_topo — Function

makevolctopo(Grid::LaMEM_grid; center::Array{Float64, 1}, height::Float64, radius::Float64, crater::Float64, base=0.0m, background=nothing)

Creates a generic volcano topography (cones and truncated cones)

Parameters

Grid - LaMEM grid (created by readLaMEMinputfile)
center - x- and -coordinates of center of volcano
height - height of volcano
radius - radius of volcano

Optional Parameters

crater - this will create a truncated cone and the option defines the radius of the flat top
base - this sets the flat topography around the volcano
background - this allows loading in a topography and only adding the volcano on top (also allows stacking of several cones to get a volcano with different slopes)

Example

Cylinder with constant phase and temperature:

julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Topo = make_volc_topo(Grid, center=[0.0,0.0], height=0.4, radius=1.5, crater=0.5, base=0.1)
CartData
    size    : (33, 33, 1)
    x       ϵ [ -3.0 : 3.0]
    y       ϵ [ -2.0 : 2.0]
    z       ϵ [ 0.1 : 0.4]
    fields  : (:Topography,)
  attributes: ["note"]
julia> Topo = make_volc_topo(Grid, center=[0.0,0.0], height=0.8, radius=0.5, crater=0.0, base=0.4, background=Topo.fields.Topography)
CartData
    size    : (33, 33, 1)
    x       ϵ [ -3.0 : 3.0]
    y       ϵ [ -2.0 : 2.0]
    z       ϵ [ 0.1 : 0.8]
    fields  : (:Topography,)
  attributes: ["note"]
julia> write_paraview(Topo,"VolcanoTopo")           # Save topography to paraview
Saved file: VolcanoTopo.vts

source

GeophysicalModelGenerator.ConstantTemp — Type

ConstantTemp(T=1000)

Sets a constant temperature inside the box

Parameters

T : the value

source

GeophysicalModelGenerator.LinearTemp — Type

LinearTemp(Ttop=0, Tbot=1000)

Set a linear temperature structure from top to bottom

Parameters

Ttop : the value @ the top
Tbot : the value @ the bottom

source

GeophysicalModelGenerator.HalfspaceCoolingTemp — Type

HalfspaceCoolingTemp(Tsurface=0, Tmantle=1350, Age=60, Adiabat=0)

Sets a halfspace temperature structure in plate

Parameters

Tsurface : surface temperature [C]
Tmantle : mantle temperature [C]
Age : Thermal Age of plate [Myrs]
Adiabat : Mantle Adiabat [K/km]

source

GeophysicalModelGenerator.SpreadingRateTemp — Type

SpreadingRateTemp(Tsurface=0, Tmantle=1350, Adiabat=0, MORside="left",SpreadingVel=3, AgeRidge=0, maxAge=80)

Sets a halfspace temperature structure within the box, combined with a spreading rate (which implies that the plate age varies)

Parameters

Tsurface : surface temperature [C]
Tmantle : mantle temperature [C]
Adiabat : Mantle Adiabat [K/km]
MORside : side of the box where the MOR is located ["left","right","front","back"]
SpreadingVel : spreading velocity [cm/yr]
AgeRidge : thermal age of the ridge [Myrs]
maxAge : maximum thermal Age of plate [Myrs]

Note: the thermal age at the mid oceanic ridge is set to 1 year to avoid division by zero

source

GeophysicalModelGenerator.LithosphericTemp — Type

LithosphericTemp(Tsurface=0.0, Tpot=1350.0, dTadi=0.5,
                    ubound="const", lbound="const, utbf = 50.0e-3, ltbf = 10.0e-3,
                    age = 120.0, dtfac = 0.9, nz = 201,
                    rheology = example_CLrheology()
                )

Calculates a 1D temperature profile [C] for variable thermal parameters including radiogenic heat source and linearly interpolates the temperature profile onto the box. The thermal parameters are defined in rheology and the structure of the lithosphere is define by LithosphericPhases().

Parameters

Tsurface : surface temperature [C]
Tpot : potential mantle temperature [C]
dTadi : adiabatic gradient [K/km]
ubound : Upper thermal boundary condition ["const","flux"]
lbound : Lower thermal boundary condition ["const","flux"]
utbf : Upper thermal heat flux [W/m]; if ubound == "flux"
ltbf : Lower thermal heat flux [W/m]; if lbound == "flux"
age : age of the lithosphere [Ma]
dtfac : Diffusion stability criterion to calculate T_age
nz : Grid spacing for the 1D profile within the box
rheology : Structure containing the thermal parameters for each phase [default example_CLrheology]

source

GeophysicalModelGenerator.ConstantPhase — Type

ConstantPhase(phase=1)

Sets a constant phase inside the box

Parameters

phase : the value

source

GeophysicalModelGenerator.compute_phase — Function

Phase = compute_phase(Phase, Temp, X, Y, Z, s::LithosphericPhases, Ztop)

Phase = compute_phase(Phase, Temp, Grid::AbstractGeneralGrid, s::LithosphericPhases)

This copies the layered lithosphere onto the Phase matrix.

Parameters

Phase - Phase array
Temp - Temperature array
X - x-coordinate array (consistent with Phase and Temp)
Y - y-coordinate array (consistent with Phase and Temp)
Z - Vertical coordinate array (consistent with Phase and Temp)
s - LithosphericPhases
Ztop - Vertical coordinate of top of model box
Grid - Grid structure (usually obtained with readLaMEMinputfile)

source

GeophysicalModelGenerator.LithosphericPhases — Type

LithosphericPhases(Layers=[10 20 15], Phases=[1 2 3 4], Tlab=nothing )

This allows defining a layered lithosphere. Layering is defined from the top downwards.

Parameters

Layers : The thickness of each layer, ordered from top to bottom. The thickness of the last layer does not have to be specified.
Phases : The phases of the layers, ordered from top to bottom.
Tlab : Temperature of the lithosphere asthenosphere boundary. If specified, the phases at locations with T>Tlab are set to Phases[end].

source

GeophysicalModelGenerator.McKenzie_subducting_slab — Type

McKenzie_subducting_slab

Thermal structure by McKenzie for a subducted slab that is fully embedded in the mantle.

Parameters

Tsurface: Top T [C]
Tmantle: Bottom T [C]
Adiabat: Adiabatic gradient in K/km
v_cm_yr: Subduction velocity [cm/yr]
κ: Thermal diffusivity [m2/s]
it: Number iterations employed in the harmonic summation

source

GeophysicalModelGenerator.LinearWeightedTemperature — Type

LinearWeightedTemperature

Structure that defined a linear average temperature between two temperature fields as a function of distance

Parameters

w_min: Minimum weight
w_max: Maximum weight
crit_dist: Critical distance
dir: Direction of the averaging (:X, :Y or :Z)
F1: First temperature field
F2: Second temperature field

source