Subduction Zone Setup#
Author: Cian Wilson
Implementation#
Even in the simplified case, the geometry of a subduction zone is more complicated than any of the in built meshes provided to us by dolfinx. In kinematic slab models we need to describe the slab and the surrounding domain around it, including crustal layers and surface features. We are particularly interested in the dynamics near the mantle wedge corner in the sub-arc region so will likely want to employ variable resolutions, with refined cells in this area. Luckily, finite elements excel at describing these more complicated, variable resolution geometries using unstructured meshes.
Due to this complication, in this notebook we only tackle the first part of step 1, implementing a function that describes the slab surface using a spline. Since the geometry is not our primary concern here we utilize classes provided in the ../python/geometry.py module for much of the implementation.
Preamble#
Let’s start by adding the path to the modules in the python folder to the system path (so we can find the geometry module).
import sys, os
basedir = ''
if "__file__" in globals(): basedir = os.path.dirname(__file__)
sys.path.append(os.path.join(basedir, os.path.pardir, 'python'))
Let’s also load the module generated by the previous notebook to get access to the default parameters loaded there.
from sz_base import default_params, allsz_params
Then let’s load all the remaining required modules at the beginning and set up a default output directory.
import geometry as geo
from mpi4py import MPI
import matplotlib.pyplot as pl
import pathlib
if __name__ == "__main__":
output_folder = pathlib.Path(os.path.join(basedir, "output"))
output_folder.mkdir(exist_ok=True, parents=True)
Slab geometry#
In kinematic slab models the slab is typically described using a small number of points derived from seismic data which are then fitted with a spline to interpolate and extrapolate the geometry to other depths. We will use a cubic spline provided by the scipy module and wrapped for convenience in our own geometry.py python module. We need to provide the points describing the spline, some information about the resolution we desire in the mesh at various points along the spline, and information about some points that we require to be included in the spline. The most important of these are the partial and full coupling depths (partial_coupling_depth and full_coupling_depth previously loaded in default_params respectively), which will later be used as the locations where the slab becomes fully coupled to the mantle wedge. These parameters are key in determining the subduction zone thermal structure. We also include a point at slab_det_depth that we use to extract diagnostic information.
We set up the slab using the function create_slab below.
def create_slab(xs, ys, resscale, lc_depth,
**kwargs):
"""
Function to construct and return a spline object that is used to describe a subducting slab
in a kinematic-slab model of a subduction zone. Optional keyword arguments default to parameters
in the global default_params dictionary if not specified.
Arguments:
* xs - list of x points in slab spline
* ys - list of y points in slab spline (must be the same length as xs)
* resscale - resolution scale factor that multiplies all _res_fact parameters
* lc_depth - depth of lower crustal boundary ("Moho")
Keyword Arguments:
distances:
* slab_diag1_depth - starting depth of slab diagnostic region
* slab_diag2_depth - end depth of slab diagnostic region
* partial_coupling_depth - partial coupling depth on slab
* full_coupling_depth - full coupling depth on slab
* slab_det_depth - detector depth on slab
resolutions factors (that get multiplied by the resscale to get the resolutions):
* slab_diag1_res_fact - start of slab diagnostic region
* slab_diag2_res_fact - end of slab diagnostic region
* partial_coupling_depth_res_fact - partial coupling depth on slab
* full_coupling_depth_res_fact - full coupling depth on slab
surface ids:
* fault_sid - fault
* slab_sid - default slab surface id
* slab_diag_sid - diagnostic region of slab
Returns:
* slab - subduction zone slab spline instance
"""
# get input parameters
# depths
slab_diag1_depth = kwargs.get('slab_diag1_depth', default_params['slab_diag1_depth'])
slab_diag2_depth = kwargs.get('slab_diag2_depth', default_params['slab_diag2_depth'])
partial_coupling_depth = kwargs.get('partial_coupling_depth', default_params['partial_coupling_depth'])
full_coupling_depth = kwargs.get('full_coupling_depth', default_params['full_coupling_depth'])
slab_det_depth = kwargs.get('slab_det_depth', default_params['slab_det_depth'])
# resolutions
slab_diag1_res = kwargs.get('slab_diag1_res_fact', default_params['slab_diag1_res_fact'])*resscale
slab_diag2_res = kwargs.get('slab_diag2_res_fact', default_params['slab_diag2_res_fact'])*resscale
partial_coupling_depth_res = kwargs.get('partial_coupling_depth_res_fact', default_params['partial_coupling_depth_res_fact'])*resscale
full_coupling_depth_res = kwargs.get('full_coupling_depth_res_fact', default_params['full_coupling_depth_res_fact'])*resscale
# surface ids
fault_sid = kwargs.get('fault_sid', default_params['fault_sid'])
slab_sid = kwargs.get('slab_sid', default_params['slab_sid'])
slab_diag_sid = kwargs.get('slab_diag_sid', default_params['slab_diag_sid'])
# set up resolutions along the slab depending on depth
# high resolution at shallow depths, lower resolution below the "diagnostic"
# region required in the benchmark case
# FIXME: these are currently hard-coded relative to the resolutions specified at the partial and full coupling
# depths for simplicity but could be separate parameters
res = [partial_coupling_depth_res if y >= -slab_diag2_depth else 3*full_coupling_depth_res for y in ys]
# set up the surface ids for the slab depending on depth
# above the "Moho" use fault_sid
# in the diagnostic region use the slab_diag_sid
# everywhere else use the default slab_sid
sids = []
for y in ys[1:]:
if y >= -lc_depth:
sid = fault_sid
elif y >= -slab_diag1_depth:
sid = slab_sid
elif y >= -slab_diag2_depth:
sid = slab_diag_sid
else:
sid = slab_sid
sids.append(sid)
# set up the slab spline object
slab = geo.SlabSpline(xs, ys, res=res, sid=sids, name="Slab")
assert full_coupling_depth > partial_coupling_depth
# adding the coupling depths may or may not be necessary
# depending on if they were included in the slab spline data already or not
# the slab class should ignore them if they aren't necessary
slab.addpoint(partial_coupling_depth, "Slab::PartialCouplingDepth",
res=partial_coupling_depth_res,
sid=slab_diag_sid)
slab.addpoint(full_coupling_depth, "Slab::FullCouplingDepth",
res=full_coupling_depth_res,
sid=slab_diag_sid)
# add the slab detector point
slab.addpoint(slab_det_depth, "Slab::DetectorPoint",
res=full_coupling_depth_res,
sid=slab_diag_sid)
# and return it
return slab
Demonstration - Benchmark#
Describing the slab geometry only takes a few non-default parameters, which are relatively simple in the simplified benchmark geometry.
Although the resolution of our mesh is going to vary across the domain we will use a resolution scale factor resscale to scale the resolution globally, while different points in the domain retain the same default relative resolutions. So a large resscale means low resolution and a small resscale means high resolution.
Computational cost
Setting the resscale too low will result in a computationally expensive simulation, especially in the non-linear case, that may need to be run locally rather than online.
if __name__ == "__main__":
resscale = 5.0
The benchmark slab geometry is rather simple, just consisting of a straight line with 2:1 horizontal distance to depth ratio, extending to 200km depth. We can therefore just provide the spline with a series of linearly related points xs and ys.
if __name__ == "__main__":
# points in slab (just linear)
xs = [0.0, 140.0, 240.0, 400.0]
ys = [0.0, -70.0, -120.0, -200.0]
To get the surface ids on the slab correct we also have to provide the lower crustal depth lc_depth. As this is a case dependent parameter it is not provided in default_params. For the benchmark cases it is at 40km depth.
if __name__ == "__main__":
lc_depth = 40
Providing these parameters we can create our slab geometry.
if __name__ == "__main__":
slab = create_slab(xs, ys, resscale, lc_depth)
/__w/fenics-sz/fenics-sz/notebooks/../python/geometry.py:214: RuntimeWarning: The iteration is not making good progress, as measured by the
improvement from the last ten iterations.
return opt.fsolve(lambda x: self.x2delu(x, x0=x0)-delu, x0, fprime=lambda x: [self.du(x)])[0]
We can double check that it looks as expected by plotting the slab, though in the benchmark case this is not very interesting!
def plot_slab(slab):
interpx = [curve.points[0].x for curve in slab.interpcurves]+[slab.interpcurves[-1].points[1].x]
interpy = [curve.points[0].y for curve in slab.interpcurves]+[slab.interpcurves[-1].points[1].y]
fig = pl.figure()
ax = fig.gca()
ax.plot(interpx, interpy)
ax.set_xlabel('x (km)')
ax.set_ylabel('y (km)')
ax.set_aspect('equal')
ax.set_title('Slab Geometry')
return fig
pl.savefig(output_folder / 'sz_slab_benchmark.png')
if __name__ == "__main__":
fig = plot_slab(slab)
fig.gca().set_title('Benchmark Slab Geometry')
fig.savefig(output_folder / 'sz_slab_benchmark.png')
Demonstration - Alaska Peninsula#
Since the benchmark geometry is not very interesting we can also demonstrate the create_slab function on a more interesting case from the global suite, “01_Alaska_Peninsula”.
if __name__ == "__main__":
resscale = 5.0
szdict_ak = allsz_params['01_Alaska_Peninsula']
slab_ak = create_slab(szdict_ak['xs'], szdict_ak['ys'], resscale, szdict_ak['lc_depth'])
and plot the resulting geometry
if __name__ == "__main__":
fig = plot_slab(slab_ak)
fig.gca().set_title('Alaska Peninsula Slab Geometry')
pl.savefig(output_folder / 'sz_slab_ak.png')
While this won’t be fully tested until we compare against existing simulations, create_slab appears to be working and we can move on to describing the rest of the subduction zone geometry in the next notebook.
Finish up#
Convert this notebook to a python script (making sure to save first) so that we can access it in subsequent notebooks.
if __name__ == "__main__" and "__file__" not in globals():
from ipylab import JupyterFrontEnd
app = JupyterFrontEnd()
app.commands.execute('docmanager:save')
!jupyter nbconvert --NbConvertApp.export_format=script --ClearOutputPreprocessor.enabled=True sz_slab.ipynb
[NbConvertApp] Converting notebook sz_slab.ipynb to script
[NbConvertApp] Writing 11643 bytes to sz_slab.py