<jupyter_start><jupyter_text>Modeling Source-to-Sink systems using FastScape: 13. sea level variations ![Namibian margin](Margin.jpg "Namibian margin")<jupyter_code>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<jupyter_output><empty_output><jupyter_text>In this last experiment we will investigate how cyclic sea level variations are recorded in the stratigraphic record.<jupyter_code>from fastscape.models import marine_model

from fastscape.processes import (BlockUplift)

sink_model = (marine_model.
              drop_processes('diffusion').
              drop_processes('init_topography').
              drop_processes('uplift').
              update_processes({'uplift': BlockUplift}))

#sink_model.visualize(show_variables=True)

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-4
u = np.zeros((ny,nx))
ylim = 2*yl/(nx-1)
u = np.where(y<ylim, u0, 0)

h0 = np.where(y>ylim, -1000*(y-ylim)/(yl - ylim),0)
h0 = np.where(h0<-400, -400, h0)
nstep = 201 # total number of steps
neq = 101 # number of steps to reach steady-state

teq = 1e7 # time to reach steady-state
period = 1e6 # period of climatic forcing
tfinal = teq + 5*period # final time

# Here we build the time array (note that not all time steps are of the same length)
tim1 = np.linspace(0,teq,101) 
tim2 = np.linspace(teq + period/10, tfinal, 100)
tim = np.concatenate((tim1,tim2))

# build precipitation array
#sea_level = np.where(tim>teq, 100*np.sin(2*np.pi*(tim-teq)/period), 0)
sea_level = 100*np.sin(2*np.pi*(tim-teq)/period)

sea_xr = xr.DataArray(data=sea_level, dims=['time'])

fig, ax = plt.subplots(nrows = 1, ncols = 1, sharex=False, sharey=True, figsize=(12,7))

ax.plot(tim, sea_level)


# %create_setup sink_model --default --verbose
import xsimlab as xs

ds_in = xs.create_setup(
    model=sink_model,
    clocks={'time': tim,
           'strati': tim[::10]},
    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'],
        # 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-4,
        # detached bedrock transport/deposition coefficient
        'spl__g_coef_bedrock': 1,
        # soil (sediment) transport/deposition coefficient
        'spl__g_coef_soil': 1,
        # sand/silt ratio of continental sediment source
        'marine__ss_ratio_land': 0.5,
        # surface (reference) porosity of sand
        'marine__porosity_sand': 0.49,
        # surface (reference) porosity of silt
        'marine__porosity_silt': 0.63,
        # e-folding depth of exp. porosity curve for sand
        'marine__e_depth_sand': 3.5e3,
        # e-folding depth of exp. porosity curve for silt
        'marine__e_depth_silt': 2e3,
        # diffusivity (transport coefficient) for sand
        'marine__diffusivity_sand': 1e1,
        # diffusivity (transport coefficient) for silt
        'marine__diffusivity_silt': 1e1,
        # mean depth (thickness) of marine active layer
        'marine__layer_depth': 100,
        # surface topography elevation
        'topography__elevation': h0,
        # sea level (elevation)
        'sea__level': sea_xr,
        # horizon freezing (deactivation) time
        'strati__freeze_time': tim,
        # uplift rate
        'uplift__rate': u,
    },
    output_vars={'topography__elevation': 'time',
                'drainage__area': 'time',
                'strati__elevation': 'strati'}
)

with xs.monitoring.ProgressBar():
    ds_out = ds_in.xsimlab.run(model=sink_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()
fig, ax = plt.subplots(figsize=(12,8))

nout = 101
for iout in range(nout-1, -1, -1):
#    ds_out.strati__elevation.isel(strati=-1).isel(horizon=iout).sel(x=xl/2)[ds_out.y>ylim].plot()
    ds_out.strati__elevation.isel(strati=-1).isel(horizon=iout).mean('x')[ds_out.y>ylim].plot()
fig, ax = plt.subplots(figsize=(12,8))

nout = 101
for iout in range(nout-1, -1, -1):
    ds_out.strati__elevation.isel(strati=-1).isel(horizon=iout).sel(y=yl/2).plot()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12779487179487178, 'fraction_numerical': 0.028923076923076923, 'mean_word_length': 3.485740570377185, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '193668', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1873, 'n_tokens_neox': 1746, 'n_words': 390, 'path': 'notebooks/MarineSink/FastScape_S2S_marine_3.ipynb', 'repo_name': 'fastscape-lem/s2s-future-dragonstone', 'size': '4875'}