Analysis and visualization of the NA direct-search results (II)

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

Import the output data

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()

Evolution of the synthetic hill: soil thickness and denudation rates

In this section we show the evolution of the synthetic hill both in terms of space-averaged denudation rate and spaced-averaged soil thickness. Here below we define the functions to calculate these time series.


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]:
30 kyr average (mm/kyr) 10 kyr average (mm/kyr) last kyr (mm/kyr)
cluster 1 5.115053 4.497795 4.594593
cluster 2 0.942027 0.742378 0.801782
cluster 3 46.713568 28.524097 25.797723
cluster 4 47.350727 26.847515 23.900095
cluster 5 59.228521 26.788265 22.656279

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.

Soil thickness distribution at characteristic times

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.

Information about this notebook

Authors: B. Bovy, J. Braun, A. Demoulin

F.R.S-FNRS | UGPQ-Ulg | ISTerre-Université de Grenoble

Softwares and versions used:


In [11]:
%load_ext version_information
%version_information numpy, scipy, pandas, tables, matplotlib, seaborn


Out[11]:
SoftwareVersion
Python2.7.9 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
IPython2.3.1
OSLinux 3.2.0 0.bpo.4 amd64 x86_64 with debian 6.0.9
numpy1.9.1
scipy0.14.0
pandas0.15.2
tables3.1.1
matplotlib1.4.2
seaborn0.5.1
Mon Jan 12 17:37:51 2015 CET

In [11]: