<jupyter_start><jupyter_text>Modeling Source-to-Sink systems using FastScape: 9. Cyclic variations in climate/precipitation ![Lannemezan Fan](LannemezanFan.jpg "Lannemezan Fan")<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 experiment, we will focus on climatic/precipitation variations that only affect the basin, not the source. In this way we can appreciate the effect of increased/decreased precipitation on the transport capacity of the depositional/transit system.<jupyter_code>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(show_inputs=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-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))<jupyter_output><empty_output><jupyter_text>To make sure that the source area is not affected by the precipitation, we will keep $K_f$ in the source constant. For this we need to create a 3D array for $K_f$ that contains both a spatial and temporal variations.<jupyter_code>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
precip = np.where(tim>teq, 1 + 0.5*np.sin(2*np.pi*(tim-teq)/period), 1)

# build Kf array and transform it into an xarray of dimension 'time' and 'space'
Kf = 1e-5
m = 0.4
Kf_tim = np.where(tim>teq, Kf*precip**m, Kf)
Kf_tim_space = np.broadcast_to(Kf_tim,(ny,nx,len(tim))).copy()
Kf_tim_space[:,:2,:] = Kf
G_tim = np.where(tim>teq, 0.5/precip, 0.5)

# create xarrays to provide adequate dimensions for FastScape
Kf_xr = xr.DataArray(data=Kf_tim_space.transpose(), dims=['time','y', 'x'])
G_xr = xr.DataArray(data=G_tim, dims=['time'])

# plots the variations of $G$ with time
fig, ax = plt.subplots(nrows = 1, ncols = 1, sharex=False, sharey=True, figsize=(12,7))

ax.plot(tim, G_tim)

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

ds_in = xs.create_setup(
    model=transit_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'],
        'uplift__rate': u,
        # MFD partioner slope exponent
        'flow__slope_exp': 1,
        # drainage area exponent
        'spl__area_exp': m,
        # slope exponent
        'spl__slope_exp': 1,
        # bedrock channel incision coefficient
        'spl__k_coef_bedrock': Kf_xr,
        # soil (sediment) channel incision coefficient
        'spl__k_coef_soil': Kf_xr,
        # detached bedrock transport/deposition coefficient
        'spl__g_coef_bedrock': G_xr,
        # soil (sediment) transport/deposition coefficient
        'spl__g_coef_soil': G_xr,
        # surface topography elevation
        'topography__elevation': np.random.random((ny,nx)),
        # horizon freezing (deactivation) time
        'strati__freeze_time': tim,
    },
    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[99])

app.show()<jupyter_output><empty_output><jupyter_text>We see that the climatic variations are imprinted in the stratigraphy but not everywhere in the basin: the variations are largest in the vicinity of the main gorge coming out of the mountain where the principal fan forms. Away from the fan, the perturbations are not stored.Also th ebasin needs to be relatively filled for the climatic signal to be stored (i.e., the amplitude of the sugnal is much smaller in the early stages of development of the basin).<jupyter_code>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()
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=ylim*3).plot()<jupyter_output><empty_output><jupyter_text>We can also compute the flux coming out of the basin (i.e. across the model boundary) to estimate, by comparison with the precipitation signal, how the climatic signal has been "filtered" or not by the basin/depositional system.<jupyter_code>nstep = len(ds_out.time)

flux = [0]
sumtop0 = ds_out.topography__elevation.isel(time=0).where(ds_out.y>=ylim).sum()
for step in range(1,nstep):
    sumtop = ds_out.topography__elevation.isel(time=step).where(ds_out.y>=ylim).sum()
    flux.append(
       (sumtop0 - sumtop)/
        (ds_out.time.values[step] - ds_out.time.values[step-1])
           )
    sumtop0 = sumtop

total_area = ds_out.grid__shape[0].values*ds_out.grid__shape[1].values
flux0 = ds_out.uplift__rate.mean().values*total_area
flux = flux/flux0<jupyter_output><empty_output><jupyter_text>Here we plot the normalized excess/default of the out-going flux (at steady-state it should be equal to 0)<jupyter_code>fig, ax = plt.subplots(nrows = 1, ncols = 1, sharex=False, sharey=True, figsize=(12,7))

ax.plot(tim, flux, label='flux')
ax.plot(tim, precip, label='precip')
ax.legend()<jupyter_output><empty_output><jupyter_text>Finally, we can compute the gain (ratio of relative amplitudes between response and forcing) and time lag.<jupyter_code>mid = 101
amp_flux = flux[mid:].max() - flux[mid:].min()
amp_forcing = precip[mid:].max() - precip[mid:].min()

print('forcing:',amp_forcing,'response:', amp_flux)


print('time lag:',(tim[np.argmax(precip[180:])+180] - tim[np.argmax(flux[180:])+180])/period)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12463851937536148, 'fraction_numerical': 0.019809138230190863, 'mean_word_length': 4.3371913580246915, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 35, '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': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '48125', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2566, 'n_tokens_neox': 2427, 'n_words': 615, 'path': 'notebooks/BasinTransit/FastScape_S2S_Transit_4.ipynb', 'repo_name': 'fastscape-lem/s2s-future-dragonstone', 'size': '6916'}