<jupyter_start><jupyter_text>Cellular model of weathering starting from a flat surfaceThis notebook runs the Grain Hill model with an initially flat surface and "walls" on either boundary, to explore what happens when you have weathering without any net erosion or deposition.Start with some imports...<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Next, set up model parameters in a dictionary:<jupyter_code># Dictionary for parameters
params = {}

num_cols = 401
params['number_of_node_columns'] = num_cols
num_rows = 41
params['number_of_node_rows'] = num_rows
params['disturbance_rate'] = 0.01
params['uplift_interval'] = 1.0e10
params['weathering_rate'] = 0.01
params['run_duration'] = 4000.0
params['show_plots'] = False
params['plot_interval'] = 40.0
params['output_interval'] = 1.1 * params['run_duration']
params['report_interval'] = 20.0
params['settling_rate'] = 220000000.0
params['friction_coef'] = 1.0
params['rock_state_for_uplift'] = 8
params['opt_rock_collapse'] = False
params['plot_every_transition'] = False<jupyter_output><empty_output><jupyter_text>To set up the starting grid, we'll use the admittedly clunky trick of creating a Landlab hex grid simply for the sake of knowing the $(x,y)$ coordinates of the nodes in the grid we will *actually* use in the Grain Hill model. We do it this way because we need to have the node-state values already defined when we create the Grain Hill model instance. The setup here has the grid half full of rock cells, with high walls of rock on the left and right boundaries.<jupyter_code>from landlab import HexModelGrid

grid = HexModelGrid(num_rows, num_cols, orientation='vert', shape='rect')
ns = np.zeros(grid.number_of_nodes, dtype=int)
bottom_half = np.where(grid.y_of_node < 0.5 * num_rows)[0]
ns[bottom_half] = 8
left_edge = np.where(grid.x_of_node == 0.0)[0]
ns[left_edge] = 8
right_edge = np.where(grid.x_of_node > 0.866 * num_cols - 1.0)[0]
ns[right_edge] = 8
params['initial_state_grid'] = ns

# Keep a copy of the original node state grid for reuse
ns_original = ns.copy()<jupyter_output><empty_output><jupyter_text>Plot the initial condition<jupyter_code>ax = grid.hexplot(ns)
plt.axis('equal')
plt.axis('off')<jupyter_output><empty_output><jupyter_text>Because we'll want to run the model several times, using different disturbance rates, we'll create a function to do so. It takes *params* as an argument, and through this dict we can set our desired disturbance and weathering rates.<jupyter_code>def run_flat_surface_weathering(params, niter):
    """Run weathering and disturbance starting with a flat surface."""
    import grain_hill_as_class

    # instantiate a GrainHill model
    num_rows = params['number_of_node_rows']
    num_cols = params['number_of_node_columns']
    grain_hill = grain_hill_as_class.GrainHill((num_rows, num_cols), **params)
    
    dt = params['run_duration'] / float(niter)

    soil_thickness = np.zeros(niter + 1)
    
    for i in range(niter):
        
        grain_hill.run_for(dt)

        # compute and write the results
        (elev_profile, soil) = grain_hill.get_profile_and_soil_thickness(grain_hill.grid, 
                                                                 grain_hill.ca.node_state)
        #plot_hill(grain_hill.ca.grid)
        soil_thickness[i+1] = np.mean(soil)

    return np.arange(0.0, dt * niter + 2, dt), soil_thickness<jupyter_output><empty_output><jupyter_text>Now let's give it a try, initially with zero disturbance rate:<jupyter_code>params['disturbance_rate'] = 0.0
(t, H0) = run_flat_surface_weathering(params, 1000)
plt.plot(t, H0)<jupyter_output><empty_output><jupyter_text>Now, we'll do a set of experimental runs in which the ratio of disturbance to weathering rate ranges 0.1, 1.0, 10.0 (this takes a few minutes to run).<jupyter_code>params['disturbance_rate'] = 0.001
params['initial_state_grid'] = ns_original.copy()
(t, H0p001) = run_flat_surface_weathering(params, 1000)
w = 0.01
plt.plot(t * w, H0)
plt.plot(t * w, H0p001)
plt.xlabel('Dimensionless time (tw)')
plt.ylabel('Dimensionless soil thickness (H / \delta)')
params['disturbance_rate'] = 0.01
params['initial_state_grid'] = ns_original.copy()
(t, H0p01) = run_flat_surface_weathering(params, 1000)
plt.plot(t * w, H0)
plt.plot(t * w, H0p001)
plt.plot(t * w, H0p01)
plt.xlabel('Dimensionless time (tw)')
plt.ylabel('Dimensionless soil thickness (H / \delta)')
#plt.legend('d/w = 0', 'd/w = 0.1', 'd/w = 1', 'd/w = 10')
params['disturbance_rate'] = 0.1
params['initial_state_grid'] = ns_original.copy()
(t, H0p1) = run_flat_surface_weathering(params, 1000)
plt.plot(t * w, H0)
plt.plot(t * w, H0p001)
plt.plot(t * w, H0p01)
plt.plot(t * w, H0p1)
plt.xlabel('Dimensionless time (tw)')
plt.ylabel('Dimensionless soil thickness (H / \delta)')<jupyter_output><empty_output><jupyter_text>Now to compare with analytical theory for thickening over time, which predicts a log relation:$h = \ln \left[ w t + 1 \right]$<jupyter_code>plt.plot(t * w, H0)
plt.plot(t * w, H0p001)
plt.plot(t * w, H0p01)
plt.plot(t * w, H0p1)
h_pred = np.log(w*t + 1.0)
plt.plot(t * w, h_pred, 'k--')
plt.xlabel('Dimensionless time (tw)')
plt.ylabel('Dimensionless soil thickness (H / \delta)')<jupyter_output><empty_output><jupyter_text>Clearly the cellular and exponential-analytical models differ at long times, with the former predicting a strong asymptote at $< 2\delta$ and the latter predicting continued growth in the thickness of the soil mantle. Which is more correct?<jupyter_code>h_pred = np.log(w*t[:500] + 1.0)
p_analytical = plt.plot(t[:500] * w, h_pred, 'c', label='log growth')
p0p1 = plt.plot(t[:500] * w, H0p1[:500], 'k', label='d/w = 10')
p0p01 = plt.plot(t[:500] * w, H0p01[:500], 'k-.', label='d/w = 1')
p0p001 = plt.plot(t[:500] * w, H0p001[:500], 'k:', label='d/w = 0.1')
p0 = plt.plot(t[:500] * w, H0[:500], 'k--', label='d/w = 0')
plt.xlabel('$t w$', fontsize=14)
plt.ylabel('$h$', fontsize=14)
plt.legend()
plt.savefig('h_vs_t_flat_surf.pdf')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1314662273476112, 'fraction_numerical': 0.03542009884678748, 'mean_word_length': 5.144736842105263, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '465496', 'license': 'MIT', 'n_tokens_mistral': 2361, 'n_tokens_neox': 2165, 'n_words': 579, 'path': 'ModelInputsAndRunScripts/WeatheringFlatSurfaceStudy/weathering_from_flat_surface.ipynb', 'repo_name': 'gregtucker/tucker_mccoy_hobley_grain_hill_manuscript', 'size': '6070'}