Concentration Tracker#

ConcentrationTrackerForDiffusion: Track mass concentration of sediment properties based on flux from DepthDependentDiffuser or DepthDependentTaylorDiffuser#

Created on Wed May 31 11:41:20 2023

@author: LaurentRoberge

class ConcentrationTrackerForDiffusion(*args, **kwds)[source]#

Bases: Component

Track the concentration of any user-defined property.

This component tracks the concentration of any user-defined property of sediment using a mass balance approach in which the concentration \(C\) is calculated as:

\[∂CH / ∂t = [-(∂q_x C_x) / ∂x - (∂q_y C_y) / ∂y] + C_br * H_brw\]

where \(H\) is sediment depth, \(q_x\) and \(q_y\) are sediment fluxed in the x and y directions, \(C_br\) is concentration in parent bedrock, and \(H_brw\) is the height of bedrock weathered into soil.

Note

This component requires a soil flux field calculated by a hillslope diffusion component and must be run after every diffusion step. Currently, this component WILL ONLY WORK IF COUPLED with the DepthDependentDiffuser or the DepthDependentTaylorDiffuser (without the dynamic timestep option).

In-situ production and decay of the material property are handled by the ConcentrationTrackerProductionDecay component.

Examples

A 1-D hillslope:

>>> import numpy as np
>>> from landlab import NodeStatus, RasterModelGrid
>>> from landlab.components import DepthDependentDiffuser
>>> from landlab.components import ConcentrationTrackerForDiffusion

>>> mg = RasterModelGrid((3, 5), xy_spacing=2.0)

>>> mg.set_status_at_node_on_edges(
...     right=NodeStatus.CLOSED,
...     top=NodeStatus.CLOSED,
...     left=NodeStatus.CLOSED,
...     bottom=NodeStatus.CLOSED,
... )
>>> mg.status_at_node[5] = NodeStatus.FIXED_VALUE
>>> mg.status_at_node.reshape(mg.shape)
array([[4, 4, 4, 4, 4],
       [1, 0, 0, 0, 4],
       [4, 4, 4, 4, 4]], dtype=uint8)

>>> mg.at_node["sediment_property__concentration"] = [
...     [0.0, 0.0, 0.0, 0.0, 0.0],
...     [0.0, 0.0, 0.0, 1.0, 0.0],
...     [0.0, 0.0, 0.0, 0.0, 0.0],
... ]
>>> mg.at_node["soil__depth"] = mg.node_x.copy()
>>> mg.at_node["bedrock__elevation"] = mg.node_x.copy()
>>> mg.at_node["topographic__elevation"] = (
...     mg.at_node["soil__depth"] + mg.at_node["bedrock__elevation"]
... )
>>> _ = mg.add_zeros("soil_production__rate", at="node")

>>> ddd = DepthDependentDiffuser(mg)
>>> ct = ConcentrationTrackerForDiffusion(mg)
>>> ddd.run_one_step(1.0)
>>> ct.run_one_step(1.0)

>>> mg.at_node["topographic__elevation"].reshape(mg.shape)
array([[ 0. ,   4.        ,   8.        ,  12.        , 16. ],
       [ 0. ,   4.11701964,   8.01583689,  11.00247875, 16. ],
       [ 0. ,   4.        ,   8.        ,  12.        , 16. ]])
>>> mg.at_node["sediment_property__concentration"].reshape(mg.shape)
array([[ 0. , 0. , 0.        , 0. , 0. ],
       [ 0. , 0. , 0.24839685, 1. , 0. ],
       [ 0. , 0. , 0.        , 0. , 0. ]])

Now, a 2-D pyramid-shaped hillslope.

>>> mg = RasterModelGrid((5, 5), xy_spacing=2.0)

>>> c = mg.add_zeros("sediment_property__concentration", at="node")
>>> c[12] = 1.0
>>> h = mg.add_full("soil__depth", 2.0, at="node")
>>> z_br = mg.add_field(
...     "bedrock__elevation",
...     8.0 - abs(4.0 - mg.node_x) - abs(4.0 - mg.node_y),
...     at="node",
... )
>>> z = mg.add_field("topographic__elevation", z_br + h, at="node")
>>> _ = mg.add_zeros("soil_production__rate", at="node")

>>> ddd = DepthDependentDiffuser(mg)
>>> ct = ConcentrationTrackerForDiffusion(mg)
>>> ddd.run_one_step(1.0)
>>> ct.run_one_step(1.0)

>>> mg.at_node["topographic__elevation"][mg.core_nodes].reshape((3, 3))
array([[ 6. , 7.13533528, 6. ],
       [ 7.13533528, 8.27067057, 7.13533528],
       [ 6. , 7.13533528, 6. ]])
>>> mg.at_node["sediment_property__concentration"][mg.core_nodes].reshape((3, 3))
array([[ 0. , 0.38079708, 0. ],
       [ 0.38079708, 1. , 0.38079708],
       [ 0. , 0.38079708, 0. ]])

And running one more step.

>>> ddd.run_one_step(1.0)
>>> ct.run_one_step(1.0)

>>> mg.at_node["topographic__elevation"][mg.core_nodes].reshape((3, 3))
array([[ 5.52060315, 6.62473963, 5.52060315],
       [ 6.62473963, 8.00144598, 6.62473963],
       [ 5.52060315, 6.62473963, 5.52060315]])
>>> mg.at_node["sediment_property__concentration"][mg.core_nodes].reshape((3, 3))
array([[ 0.09648071, 0.44750673, 0.09648071],
       [ 0.44750673, 1.        , 0.44750673],
       [ 0.09648071, 0.44750673, 0.09648071]])

Finally, the same 2D hillslope now using the DepthDependentTaylorDiffuser. Note that the timestep must be smaller than 1 to maintain stability in the diffusion calculation. Typically, one could use the dynamic timestepping option. However, here it will provide incorrect soil flux values to the ConcentrationTrackerForDiffusion, which cannot do sub-timestep calculations. Use the if_unstable=”warn” flag when instantiating the Taylor diffuser and pick a timestep that is stable.

>>> from landlab.components import DepthDependentTaylorDiffuser
>>> mg = RasterModelGrid((5, 5), xy_spacing=2.0)

>>> c = mg.add_zeros("sediment_property__concentration", at="node")
>>> c[12] = 1.0
>>> h = mg.add_full("soil__depth", 2.0, at="node")
>>> z_br = mg.add_field(
...     "bedrock__elevation",
...     8.0 - abs(4.0 - mg.node_x) - abs(4.0 - mg.node_y),
...     at="node",
... )
>>> z = mg.add_field("topographic__elevation", z_br + h, at="node")
>>> _ = mg.add_zeros("soil_production__rate", at="node")

>>> ddtd = DepthDependentTaylorDiffuser(mg, if_unstable="warn")
>>> ct = ConcentrationTrackerForDiffusion(mg)
>>> ddtd.run_one_step(0.4)
>>> ct.run_one_step(0.4)

>>> mg.at_node["topographic__elevation"][mg.core_nodes].reshape((3, 3))
array([[ 6.        , 7.30826823, 6.        ],
       [ 7.30826823, 8.61653645, 7.30826823],
       [ 6.        , 7.30826823, 6.        ]])
>>> mg.at_node["sediment_property__concentration"][mg.core_nodes].reshape((3, 3))
array([[ 0.        , 0.26436925, 0.        ],
       [ 0.26436925, 1.        , 0.26436925],
       [ 0.        , 0.26436925, 0.        ]])

References

Required Software Citation(s) Specific to this Component

CITATION

Parameters:

grid (ModelGrid) – Landlab ModelGrid object
concentration_initial (float, array, or str, optional) – Initial concentration in soil/sediment as either a scalar, array, or the name of an existing node field, -/m^3.
concentration_in_bedrock (float, array, or str, optional) – Concentration in bedrock as either a scalar, array, or the name of an existing node field, -/m^3.
concentration_from_weathering (float or array, optional) – Concentration generated during the weathering process, -/m^3. Defaults to None, which causes all weathered bedrock to retain its original parent material concentration (concentration_in_bedrock) as it weathers to soil. Use this parameter to differentiate between the concentration in weathered material compared to its parent bedrock.

__init__(grid, concentration_initial=0, concentration_in_bedrock=0, concentration_from_weathering=None)[source]#

Parameters:

grid (ModelGrid) – Landlab ModelGrid object
concentration_initial (float, array, or str, optional) – Initial concentration in soil/sediment as either a scalar, array, or the name of an existing node field, -/m^3.
concentration_in_bedrock (float, array, or str, optional) – Concentration in bedrock as either a scalar, array, or the name of an existing node field, -/m^3.
concentration_from_weathering (float or array, optional) – Concentration generated during the weathering process, -/m^3. Defaults to None, which causes all weathered bedrock to retain its original parent material concentration (concentration_in_bedrock) as it weathers to soil. Use this parameter to differentiate between the concentration in weathered material compared to its parent bedrock.

calc_concentration(dt)[source]#

Calculate change in concentration for a time period ‘dt’.

Parameters:: dt (float (time)) – The imposed timestep.

property conc_br#: Concentration in bedrock (kg/m^3).

property conc_init#: Initial concentration in soil/sediment (kg/m^3).

property conc_w#: Concentration from the weathering process (kg/m^3).

run_one_step(dt)[source]#

Run for a time step.

Parameters:: dt (float (time)) – The imposed timestep.