Source code for landlab.components.tectonics.listric_kinematic_extender

```#!/usr/bin/env python3
"""
Apply tectonic extension kinematically.

Landlab component that simulates development of an asymmetric rift on a listric
fault plane.

See notebook tutorial for theory and examples.

@author: gtucker
"""

import numpy as np

from landlab import Component, HexModelGrid, RasterModelGrid

[docs]
def dist_to_line(Px, Py, x0, y0, alpha):
"""Calculate and return the distance of point(x) (Px, Py) to the
line described by x = x0 + t cos alpha, y = y0 + t sin alpha.

Parameters
----------
Px : float
x-coordinate of point(s)
Py : float
y-coordinate of point(s)
x0 : float
x intercept of line
y0 : float
y intercept of line
alpha : float, degrees
angle of line, counter-clockwise from positive x-axis

Examples
--------
>>> np.round(dist_to_line(1, 1, 0, 0, 90), 6)
1.0
>>> np.round(dist_to_line(0, 1, 1, 0, 90), 6)
-1.0
>>> np.round(dist_to_line(1, 1, 0, 0, 0), 6)
-1.0
>>> np.round(dist_to_line(2.0**0.5, 0, 0, 0, 45), 6)
1.0
>>> np.round(dist_to_line(0, 2.0**0.5, 0, 0, 45), 6)
-1.0
"""
return np.sin(alpha_r) * (Px - x0) - np.cos(alpha_r) * (Py - y0)

[docs]
class ListricKinematicExtender(Component):
"""Apply tectonic extension kinematically to a raster or
hex grid.

The caller specifies the strike, dip, and location of the zero-surface
fault trace (i.e., where the fault plane would intersect zero elevation),
and either the (x, y) components of uniform extension velocity field,
or a link-based velocity field. The run_one_step() method calculates
advection of an output field called "hangingwall__thickness". The initial
hanginwall thickness is defined as the difference between the starting
topography field (a required input field) and a listric fault plane
that is represented mathematically as an "upside-down" saturating
exponential function that asymptotes to a caller-specified detachment
depth, representing a decollement.

Examples
--------
>>> from landlab import RasterModelGrid
>>> from landlab.components import ListricKinematicExtender
>>> grid = RasterModelGrid((3, 130), xy_spacing=10.0)
>>> lke = ListricKinematicExtender(grid, fault_x0=100.0, fault_strike=90.0)
>>> for _ in range(250):
...     lke.run_one_step(dt=2000.0)
...
>>> round(grid.at_node["hangingwall__thickness"][240])
830
"""

_name = "ListricKinematicExtender"

_time_units = "y"

_unit_agnostic = True

_info = {
"dtype": float,
"intent": "in",
"optional": True,
"units": "m/y",
},
"fault_plane__elevation": {
"dtype": "float",
"intent": "out",
"optional": False,
"units": "m",
"mapping": "node",
"doc": "Elevation of fault plane",
},
"hangingwall__thickness": {
"dtype": "float",
"intent": "out",
"optional": False,
"units": "m",
"mapping": "node",
"doc": "Thickness of material in hangingwall block",
},
"topographic__elevation": {
"dtype": "float",
"intent": "inout",
"optional": False,
"units": "m",
"mapping": "node",
"doc": "Land surface topographic elevation",
},
}

[docs]
def __init__(
self,
grid,
extension_rate_x=0.001,
extension_rate_y=0.0,
fault_dip=60.0,
fault_x0=0.0,
fault_y0=0.0,
fault_strike=45.0,
detachment_depth=1.0e4,
):
"""Deform vertically and horizontally to represent tectonic extension.

Parameters
----------
grid: RasterModelGrid
A landlab grid.
extension_rate_x: float, optional
Rate of x-directed horizontal motion of hangingwall relative to footwall
(m / y), default 0.001 m/y.
extension_rate_y: float, optional
Rate of y-directed horizontal motion of hangingwall relative to footwall
(m / y), default 0.
fault_x0: float, optional
x intercept of zero-surface fault trace, m (default 0).
fault_y0: float, optional
y intercept of zero-surface fault trace, m (default 0).
fault_strike: float, optional
Strike of zero-surface fault trace, degrees (default 45).
detachment_depth: float, optional
Depth to horizontal detachment (m), default 10 km.
List of names of fields, in addition to 'hangingwall__thickness'
Indicates whether the directions of advection are expected to remain
steady throughout a run. If True, some computation time is saved
by calculating upwind links only once.
"""
if not (isinstance(grid, RasterModelGrid) or isinstance(grid, HexModelGrid)):
raise (TypeError, "grid must be a RasterModelGrid or HexModelGrid")

super().__init__(grid)
self.initialize_output_fields()

self._elev = grid.at_node["topographic__elevation"]
self._fault_plane_elev = grid.at_node["fault_plane__elevation"]
self._hw_thick = grid.at_node["hangingwall__thickness"]

self.update_fault_plane_elevation_and_hangingwall_thickness(
grid, fault_x0, fault_y0, fault_strike, fault_dip, detachment_depth
)
grid,
extension_rate_x,
extension_rate_y,
)

[docs]
def update_fault_plane_elevation_and_hangingwall_thickness(
self, grid, fault_x0, fault_y0, fault_strike, fault_dip, detachment_depth
):
"""Initialize fields fault_plane__elevation and hangingwall__thickness.

Calculate and store the fault plane elevation at grid nodes using an
exponential function of (signed) distance to fault, with topographic
elevation as the minimum. Calculate the thickness of the hangingwall
block at grid nodes by subtracting fault plane elevation from
topographic elevation.

Parameters
----------
fault_x0 : float
x-intercept of zero-surface fault trace, m
fault_y0 : float
y-intercept of zero-surface fault trace, m
fault_strike : float
strike angle of fault trace, degrees ccw from +x
fault_dip : float
dip angle of fault at the zero elevation point, degrees
detachment_depth : float
depth to the point where the detachment is horizontal, m

Examples
--------
>>> from landlab import RasterModelGrid
>>> from landlab.components import ListricKinematicExtender
>>> grid = RasterModelGrid((3, 3), xy_spacing=1000.0)
>>> extender = ListricKinematicExtender(grid, fault_strike=90.0)
>>> round(grid.at_node["fault_plane__elevation"][4])
-1590
>>> round(grid.at_node["hangingwall__thickness"][4])
1590
"""
dist_to_fault = dist_to_line(
grid.x_of_node, grid.y_of_node, fault_x0, fault_y0, fault_strike
)
self._fault_plane_elev[:] = np.minimum(
-detachment_depth
* (1.0 - np.exp(-dist_to_fault * fault_grad / detachment_depth)),
self._elev,
)
self._hw_thick[:] = self._elev - self._fault_plane_elev

self,
grid,
extension_rate_x,
extension_rate_y,
):

non-zero values, these values are used for the advection field.
Otherwise, the field is created (if needed) and initialized by mapping
the vector components extension_rate_x and extension_rate_y onto the

List of names of fields to advect. Can be an empty list.
"hangingwall__thickness" will be added to the list if it
extension_rate_x: float
Rate of x-directed horizontal motion of hangingwall relative to footwall
(m / y).
extension_rate_y: float
Rate of y-directed horizontal motion of hangingwall relative to footwall
(m / y).
Indicates whether the directions of advection are expected to remain
steady throughout a run. If True, some computation time is saved
by calculating upwind links only once.
"""
if np.amax(np.abs(self._advec_velocity)) == 0.0:  # if no nonzero values
)
grid,
)

[docs]
def run_one_step(self, dt):
"""Apply extensional motion to grid for one time step.

Parameters
----------
dt : float
Time-step duration, y
"""