# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python [conda env:s2s-future-dragonstone]
#     language: python
#     name: conda-env-s2s-future-dragonstone-py
# ---

# # Modeling Source-to-Sink systems using FastScape: 8. Storing stratigraphy in a subsiding basin

# ![Lannemezan Fan](LannemezanFan.jpg "Lannemezan Fan")

import xsimlab as xs
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
#plt.style.use('dark_background')
# %load_ext xsimlab.ipython
import hvplot.xarray

# ## The marine model
#
# The marine model includes a diffusion model to reproduce the transport of sediment in the marine environment. It also includes a stratigraphic model. In this example, we only wish to have access to the stratigraphic model, so we are going to use the marine model but after stripping it of its marine component, the diffusion, the initial topography, and replace the uplift function by a simpler process (BlockUplift).

# +
from fastscape.models import marine_model

from fastscape.processes import (BlockUplift)

transit_model = (marine_model.
              drop_processes('diffusion').
              drop_processes('init_topography').
              drop_processes('uplift').
              drop_processes('marine').
              drop_processes('sea').
              update_processes({'uplift': BlockUplift}))

transit_model.visualize()

# -

# We introduce an uplift function that contains (in addition to the uplift in the source area) a subsidence term to create a basin. The subsidence is supposed to decrease linearly from the source/mountain edge to the base level (the opposite boundary).

# +
xl = 100e3
yl = 100e3
nx = 101
ny = 101
X = np.linspace(0,xl,nx)
Y = np.linspace(0,yl,ny)
x,y = np.meshgrid(X, Y)

u0 = 3e-2
u1 = -1e-4
u = np.zeros((ny,nx))
ylim = 2*yl/(nx-1)
u = np.where(y<ylim, u0, u1*(yl-y)/(yl-ylim))

# -

# To use the stratigraphy process, one need to specify the ages of the stratigraphic markers (in model time) in the *freeze_time* variable. We also define a new clock at which the stratigraphic information will be stored (for display purposes).

# +
# # %create_setup transit_model --default --verbose
import xsimlab as xs

ds_in = xs.create_setup(
    model=transit_model,
    clocks={'time': np.linspace(0,1e7,101),
            'strati': np.linspace(0, 1e7, 11)},
    master_clock='time',
    input_vars={
        # nb. of grid nodes in (y, x)
        'grid__shape': [ny,nx],
        # total grid length in (y, x)
        'grid__length': [yl,xl],
        # node status at borders
        'boundary__status': ['looped','looped','fixed_value','core'],
        'uplift__rate': u,
        # MFD partioner slope exponent
        'flow__slope_exp': 1,
        # drainage area exponent
        'spl__area_exp': 0.4,
        # slope exponent
        'spl__slope_exp': 1,
        # bedrock channel incision coefficient
        'spl__k_coef_bedrock': 1e-5,
        # soil (sediment) channel incision coefficient
        'spl__k_coef_soil': 1e-5,
        # detached bedrock transport/deposition coefficient
        'spl__g_coef_bedrock': 1,
        # soil (sediment) transport/deposition coefficient
        'spl__g_coef_soil': 1,
        # surface topography elevation
        'topography__elevation': np.random.random((ny,nx)),
        # horizon freezing (deactivation) time
        'strati__freeze_time': np.linspace(0,1e7,11),
    },
    output_vars={'topography__elevation': 'time',
                'drainage__area': 'time',
                'strati__elevation': 'strati'}
)

# -

with xs.monitoring.ProgressBar():
    ds_out = ds_in.xsimlab.run(model=transit_model)

# +
from ipyfastscape import TopoViz3d

app = TopoViz3d(ds_out, canvas_height=600, time_dim="time")

app.components['background_color'].set_color('lightgray')
app.components['vertical_exaggeration'].set_factor(5)
app.components['timestepper'].go_to_time(ds_out.time[-1])

app.show()
# -

# To visualize the stratigraphy, we can plot the various horizons at the last time step. First along, then across the basin.

# +
fig, ax = plt.subplots(figsize=(12,8))

nout = len(ds_out.horizon)
for iout in range(nout):
    ds_out.strati__elevation.isel(strati=-1).isel(horizon=iout).sel(x=xl/2)[ds_out.y>ylim].plot()

# +
fig, ax = plt.subplots(figsize=(12,8))

nout = len(ds_out.horizon)
for iout in range(nout):
    ds_out.strati__elevation.isel(strati=-1).isel(horizon=iout).sel(y=ylim).plot()
# -

# Anaylse this result by producing different cross sections in the $x$- and $y$-directions.



StarcoderdataJupyterScriptsDedupFiltered