Authors: B. Bovy, J. Braun, A. Demoulin.
This notebook follows a previous one in which we analysed of the outputs of the NA direct-search. We identified five clusters of models that each provide predictions that are consistent with the observed soil thickness distribution. Here, we analyse and visualize the results of forward-modelling from each cluster centroids so that we can see how these models behave differently from each other.
We start to import some Python packages and modules that will be used here.
In [1]:
import os
from collections import OrderedDict
import numpy as np
import pandas as pd
import tables
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.mlab import griddata
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns
%matplotlib inline
The directory of the model runs:
In [2]:
sim_dir = "/home/bovy/cliche/sims/hillArdn_inversion01"
We specify the output files (HDF5 format) of the five model runs...
In [3]:
hdf5_dir = os.path.join(sim_dir, "hdf5")
h5_files = OrderedDict()
for i in range(1, 6):
h5_files['cluster {}'.format(i)] = os.path.join(
hdf5_dir, "sim_cluster{}.hdf5".format(i)
)
.. and provide some functions to retreive the model run output data.
In [4]:
def get_mesh(h5_file):
"""
Return mesh arrays (x, y, ntype).
"""
with tables.openFile(h5_file) as f:
x = f.root.irrmesh.x.read()
y = f.root.irrmesh.y.read()
ntype = f.root.irrmesh.ntype.read()
return x, y, ntype
def get_zh_at_time(h5_file, time):
"""
Return z (elevation) and h (soil thickness)
arrays at the given simulation `time` (in years)
or at the end of the simulation (`time` = -1).
"""
with tables.openFile(h5_file) as f:
if time == -1:
z = f.root.variables.z.read()
h = f.root.variables.h.read()
return z, h
dt_out = f.root.params.t_params.read()[-2]
inc = int(time / dt_out)
h5_z_node = '/time_series/data_z_{}'.format(inc)
h5_h_node = '/time_series/data_h_{}'.format(inc)
z = f.get_node(h5_z_node).read()
h = f.get_node(h5_h_node).read()
return z, h
# Cell surface array is unfortuneately not saved in out hdf5.
# Get it from the mesh file
mesh_file = os.path.join(hdf5_dir, "irrmesh",
"rect_hillslope_Ardn_inv20.hdf5")
with tables.openFile(mesh_file) as f:
surface = f.root.irrmesh.surface.read()
In [5]:
def mean_denudation_time_series(h5_file, t_lim=(0., 1.2e5)):
"""
Return time series of mean denudation rate.
"""
with tables.openFile(h5_file) as f:
ntype = f.root.irrmesh.ntype.read()
mask = ntype < 1
dt_out = f.root.params.t_params.read()[-2]
inc_lim = (np.array(t_lim) / dt_out).astype('int')
inc_range = np.arange(*inc_lim)
z_h5_nodes = ['/time_series/data_z_{}'.format(i)
for i in inc_range]
z_arrays = [f.get_node(node).read()[mask]
for node in z_h5_nodes]
z = np.array(z_arrays)
diff_volumes = -1. * (z[1:] - z[:-1]) * surface[mask]
mean_denudation = np.sum(diff_volumes, axis=1)
mean_denudation /= np.sum(surface[mask])
t_range = (inc_range * dt_out)[1:]
return t_range, mean_denudation
def mean_thickness_time_series(h5_file, t_lim=(0., 1.2e5)):
"""
Return time series of mean denudation rate.
"""
with tables.openFile(h5_file) as f:
ntype = f.root.irrmesh.ntype.read()
mask = ntype < 1
dt_out = f.root.params.t_params.read()[-2]
inc_lim = (np.array(t_lim) / dt_out).astype('int')
inc_range = np.arange(*inc_lim)
h_h5_nodes = ['/time_series/data_h_{}'.format(i)
for i in inc_range]
h_arrays = [f.get_node(node).read()[mask]
for node in h_h5_nodes]
h = np.array(h_arrays)
h_mean = np.sum(h * surface[mask], axis=1)
h_mean /= np.sum(surface[mask])
t_range = (inc_range * dt_out)
return t_range, h_mean
The following code cell draws these time series over the whole simulation time-span. The color codes are the same than those used for clustering.
In [6]:
fig, axes = plt.subplots(1, 2, figsize=(12, 6))
ax1, ax2 = axes
for cluster in h5_files.keys():
t, d = mean_denudation_time_series(h5_files[cluster])
ax1.plot(t / 1e3, d * 1e3)
t, h = mean_thickness_time_series(h5_files[cluster])
ax2.plot(t / 1e3, h, label=cluster)
ax2.legend()
plt.setp(ax1, xlabel='simulation time (kyr)',
ylabel='mean denudation rate (mm/kyr)')
plt.setp(ax2, xlabel='simulation time (kyr)',
ylabel='mean soil thickness (m)')
fig.tight_layout()
We see in this figure that the model runs that we performed from cluster centroids can be further distinguished into two categories: (1) clusters 1 and 2 both present very low denudation rates and a slow evolution of space-averaged soil thickness that doesn't seems to be much responsive to changes in climatic inputs, while (2) clusters 3, 4, 5 all present much higher denudation rates and more contrasted evolution of soil thickness, which are both very responsive to climate change, and where we can see well-defined, high peaks of denudation during the transitions between cold/warm climatic stages. This cluster grouping is consistent given the higher-level clusters previously obtained from hierarchical clustering (see the dendrogram in the first notebook).
It is also worth noting that, despite the contrasted values of initial soil thickness set for each simulation, all curves of spaced-averaged soil thickness evolution shown in the figure converge towards values in the small range 0.6-0.8 m. This hightlights the low impact of the initial conditions on the state of the hillslope system at the end of the simulation. It is, however, likely that the initial soil thickness has an non-negligible influence on the values of the model parameters included in the NA search, although no simple relationship is clearly visible on the NA output ensemble (see the scatterplot matrix in the first notebook).
Now let's see how our obtained denudation rates compare with independent denudation data that have been collected by various authors in the same study-area, i.e., the main catchments in the Ardenne, and also catchements located within the surrounding Renish Massif (Schaller et al., 2001, 2002, 2004; Meyer et al., 2008). These authors inferred time and space-averaged denudation rates either from cosmogenic radio-nucleide (CRN) concentrations - measured in basal fluvial terrace or in active bedload samples - or from river load gauging. The CRN-derived rates, integrated through the last 10-30 kyr, vary between 30 and 80 mm/kyr, while the gauging-derived rates, which integrate the record over shorter periods of only a few years to decades, are about one-third of the former rates and range from 9 to 25 mm/kyr.
Comparing these denudation rates with our results is not straightforward as the former may also include the contribution of river incision, though the observations made so far tend to minimize that contribution in our case. River terrace datations made by Rixhon et. al. (2010) indeed suggest that the incision of the rivers in the Ardenne Massif occured mainly as the slow propagation of an erosion wave over time scales much longer than the 10-30 kyr considered here. Schaller et al. (2001, 2011 pers. comm.) have also observed no significant downstream increase of CRN concentrations from multiple samples taken along the streams.
Here below we calculate the time-space-averaged denudation rates from the model outputs, considering the time-spans relative to the measured rates.
In [7]:
char_den_rates = OrderedDict()
col_names = ['30 kyr average (mm/kyr)',
'10 kyr average (mm/kyr)',
'last kyr (mm/kyr)']
for cluster in h5_files.keys():
t, d = mean_denudation_time_series(h5_files[cluster])
char_den_rates[cluster] = [d[t >= 1.2e5-3e4].mean(),
d[t >= 1.2e5-1e4].mean(),
d[-1]]
char_den_rates_df = pd.DataFrame.from_dict(char_den_rates, orient='index')
char_den_rates_df *= 1e3
char_den_rates_df.columns = col_names
char_den_rates_df
Out[7]:
This table shows that the models from clusters 1 and 2 predict denudation rates much lower than the measured rates. The rates predicted by clusters 3, 4, 5 are more consistent with the measured rates, they will be further discussed.
The figure above provides useful information for characterizing the evolution of the modelled hillslope system. However, it doesn't tell us anything about the spatial configuration of the system, which might be a key in our understanding on how the behavior of the model explained the convergence of the NA search.
Here below we plot the soil distribution on the synthetic hill at two characteristic times of the simulation for each cluster. The figure above shows, at least for clusters 1-2-3, that the system dramatically responded to the climate transitions and that the erosion regime was more efficient during the cold stage than during the warm stage. It would be therefore interesting to characterize the soil thickness distribution both as (1) resulting from the erosion regime during the cold stage + cold/warm transition and (2) resulting from the response to the last cold/warm climate transition under the temperate climatic conditions. We thus here compare soil thickness distributions after 108 kyr of model run (end of cold/warm transition) and after 120 kyr (present-day).
We first define some functions for gridding the output data at the same (low) resolution than the observational data. This step is not really neccessary here but it greatly reduces the computation cost for drawing the synthetic hill (the matplotlib package is not designed for efficient 3D-plotting of high resolution data).
In [8]:
def set_grid_coord(c, grid_spacing=20.):
"""
Given any coordinate array `c`, return a gridded
coordinate array with `grid_spacing` resolution and
which is completely included in the domain of `c`.
"""
grid_coord = np.linspace(c.min(), c.max(),
c.ptp() / grid_spacing)
return grid_coord[1:-1]
def get_zh_at_time_grid(h5_file, time, grid_spacing=20.):
"""
Get z and h arrays at `time`, re-interpolated
onto a regular grid. Returns (x, y, z, h) 2-d arrays.
"""
x, y, ntype = get_mesh(h5_file)
z, h = get_zh_at_time(h5_file, time)
# grid coordinates
mask = ntype < 1
x_m, y_m = x[mask], y[mask]
x_g, y_g = set_grid_coord(x_m), set_grid_coord(y_m)
x_g2, y_g2 = np.meshgrid(x_g, y_g)
# interpolation
z_g2 = griddata(x_m, y_m, z[mask], x_g, y_g,
interp='linear')
h_g2 = griddata(x_m, y_m, h[mask], x_g, y_g,
interp='linear')
return x_g2, y_g2, z_g2.filled(), h_g2.filled()
We also provide functions to categorize the simulated soil thicknesses according to the categories defined in the soil map. This will help us to better compare the observed/modelled distributions.
In [9]:
h_categories = {
0: (0., 0.01),
1: (0.01, 0.2),
2: (0.2, 0.4),
3: (0.4, 0.8),
4: (0.8, np.Inf)
}
def get_h_mask(h, bounds):
"""
Return a mask of a soil thickness given
a soil thickness n-d array `h` and an interval
`bounds` as a (lower, upper) tuple.
The mask verify `lower < h <= upper`.
"""
lower, upper = bounds
return np.logical_and(h > lower,
h <= upper)
def get_h_categories(h):
"""
Return a categorized (integer) version
of the soil thickness n-d array `h`.
"""
h_cat = np.zeros_like(h, dtype='int')
for c, bounds in h_categories.items():
cat_mask = get_h_mask(h, bounds)
h_cat[cat_mask] = c
return h_cat
The code cell below plots the re-gridded and categorized soil thickness distributions after 108 kyr and 120 kyr for each cluster centroid. Note that the color bar is not proportional to category ranges.
In [15]:
# characteristic model run times at which the soil thickness
# distribution is plotted
run_times = (1.08e5, 1.2e5)
# define a discrete colormap for soil thickness categories
h_cat_cm = sns.cubehelix_palette(start=.5, light=1., as_cmap=True)
h_cat_lbounds = [v[0] for v in h_categories.values()] + [10.]
h_cat_norm = mpl.colors.BoundaryNorm(sorted(h_cat_lbounds),
h_cat_cm.N)
# cluster color codes
cat_clr = [mpl.colors.rgb2hex(rgb)
for rgb in sns.color_palette(n_colors=5)]
fig = plt.figure(figsize=(14, 22))
iax = 1
for c, cluster in enumerate(h5_files.keys()):
for t in run_times:
# define the subplot axe
ax = fig.add_subplot(
len(h5_files), len(run_times), iax,
projection='3d', axis_bgcolor='w'
)
# re-grid model outputs and categorize soil thickness
x_g2, y_g2, z_g2, h_g2 = get_zh_at_time_grid(
h5_files[cluster], t
)
h_cat_g2 = get_h_categories(h_g2)
clr_array = h_cat_g2 * 255 / (len(h_categories) - 1)
# 3D plot of the synthetic hill with categorized soil thickness
ax.plot_surface(x_g2, y_g2, z_g2, rstride=1, cstride=1,
antialiased=True,
facecolors=h_cat_cm(clr_array))
ax.view_init(azim=40.)
ax.set_zlim(0, 120)
# subplot tweaks
title = "{0} ({1} kyr)".format(cluster, t / 1e3)
ax.set_title(title, color=cat_clr[c])
plt.setp(ax, zlabel='elevation (m)',
xlabel='x (m)', ylabel='y (m)')
iax += 1
# colorbar and subplot adjustements
fig.subplots_adjust(bottom=0.1)
h_cat_m = plt.cm.ScalarMappable(cmap=h_cat_cm,
norm=h_cat_norm)
h_cat_m.set_array(h_g2)
cax = fig.add_axes([0, 0.05, 1, 0.1])
cax.axis('off')
clb = fig.colorbar(h_cat_m, ax=cax,
extend='max',
orientation='horizontal')
clb.set_label("categorized soil thickness (m)")
We see again in this figure that clusters 1-2 and clusters 3-4-5 form two categories, here with each very similar soil thickness distributions. The distributions haven't much evovled from 108 kyr to 120 kyr for clusters 1-2, while it changed more drastically for clusters 3-4-5.
For this last group of clusters, we can qualify the hillslope system at 108 kyr as production-limited on a large part of the synthetic hill covered by shallow soils, while the system at 120 kyr is rather transport-limited with a globally thicker soil cover.
Authors: B. Bovy, J. Braun, A. Demoulin
F.R.S-FNRS | UGPQ-Ulg | ISTerre-Université de Grenoble
This work is licensed under a Creative Commons Attribution 4.0 International License.
Softwares and versions used:
In [11]:
%load_ext version_information
%version_information numpy, scipy, pandas, tables, matplotlib, seaborn
Out[11]:
In [11]: