# Kriging Alluvial Thicknesses in Valley Bottoms Using Nonstationary Geometric Anisotropies

Jean-Louis Grimaud, Nicolas Desassis, Diana Chourio-Camacho, Didier Renard, Mike Pereira, Fabien Ors, Hélène Tissoux, Paul Bessin, Mark Noble

► To cite this version:

Jean-Louis Grimaud, Nicolas Desassis, Diana Chourio-Camacho, Didier Renard, Mike Pereira, et al.. Kriging Alluvial Thicknesses in Valley Bottoms Using Nonstationary Geometric Anisotropies. Mathematical Geosciences, 2025, 10.1007/s11004-025-10200-w. hal-05142486

HAL Id: hal-05142486

https://hal.science/hal-05142486v1

Submitted on 3 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Kriging Alluvial Thicknesses in Valley Bottoms Using Nonstationary Geometric Anisotropies

Jean-Louis Grimaud<fnref n="1" />, ∗ · Nicolas Desassis<fnref n="1" /> · Diana Chourio-Camacho<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> · Didier Renard<fnref n="1" /> · Mike Pereira<fnref n="1" /> · Fabien Ors<fnref n="1" /> · Hélène Tissoux<fnref n="3" />,<fnref n="4" /> · Paul Bessin<fnref n="5" /> · Mark Noble<fnref n="1" />

Received: 11 March 2024 / Accepted: 3 May 2025
© The Author(s) 2025

Abstract

Modeling the geometry and volumes of alluvial infill at the bottom of river valleys is essential for a better understanding of river dynamics and alluvial storage, which buffers sedimentary signal transmission from source to sink. This study introduces kriging methods for predicting alluvial thickness geometries in the valley bottom using random fields with nonstationary geometric anisotropies. These anisotropies are constructed based on the contours of the alluvial plain. The two methods (one based on covariance and the other on Stochastic Partial Differential Equations (SPDE)) are compared to kriging with stationary anisotropies. The main parameters are estimated by maximum likelihood and the different methods are tested using cross-validation. The methods incorporating nonstationary anisotropies exhibit lower uncertainties than the stationary anisotropy approach and produce patterns that align more realistically with river geomorphic and sedimentologic knowledge. Evidence of braided patterns, riflemen pool geometry related to meandering, and scouring at confluences is observed at the scale of the valley width. Finally, it is estimated that about 2.30 km<fnref n="3" /> of alluvium is currently stored in the studied Oise valley bottom. Although the method was applied to an alluvial setting, it could be easily adapted to other geological contexts where

<figure><img src="image_2.png" /></figure>

<footnote>Jean-Louis Grimaud
jean-louis.grimaud@minesparis.psl.eu</footnote>

<footnote n="1">1 PSL University/ MINES Paris/ Centre de Géosciences, 35 rue St Honoré, 77305 Fontainebleau
Cedex, France</footnote>

<footnote n="2">2 Sorbonne Université, CNRS, EPHE, UMR 7619 METIS, 75005 Paris, France</footnote>

<footnote n="3">³ BRGM, 3 Avenue Claude Guillemin, BP 36009, 45060 Orléans, France</footnote>

<footnote n="4">4 Histoire Naturelle de l'Homme Préhistorique (HNHP), UMR 7194, Département « Homme et Environnement », Muséum National d'Histoire Naturelle, 1 rue René Panhard, 75013 Paris, France</footnote>

<footnote n="5">5 Laboratoire de Planétologie et Géosciences, LPG UMR 6112, CNRS, Le Mans Université, Univ Angers, Nantes Université, Avenue Olivier Messiaen, 72085 Le Mans, France</footnote>
anisotropy fields can be estimated, such as faulted and folded structures, supergene mineralization or environmental context to track pollutants dissemination in rivers.

**Keywords** Alluvium · Kriging · Nonstationary anisotropies

# 1 Introduction

Alluvial valleys shaped by rivers currently support the residences and infrastructures of billions of people, provide critical ecological and economic services (Best and Darby 2020), and contain paleoclimatic archives and remnants of early human settlements (Pastre et al. 2000; Grimaud et al. 2023). The volume of alluvium stored in alluvial plains reflects the landscape's capacity to temporarily store alluvium in the transfer zone, thereby buffering sedimentary signal transmission from source to sink (Castell-tort and Driessche, 2003; Blöthe and Korup 2013). Additionally, reconstructing the geometry of river valleys and volumes of alluvium stored at their bottom is crucial for a comprehensive understanding of river dynamics and alluvial storage during the Late Quaternary (e.g., Pastre et al. 2000; Antoine et al. 2007; Bridgeland and Westaway, 2008). Recent studies have demonstrated that the geometry and incision dynamics of fluvial corridors are influenced not only by external forcings but also by bedrock lithology and total valley height (e.g., Grimaud et al. 2014; Langston and Temme 2019; Tofelde et al. 2022). A wide variety of tools are now available for measuring the planform geometry of valleys and terraces using digital elevation models (e.g., Clubb et al. 2017, 2022). In addition, complementary approaches include reconstructing the full 3D geometry of alluvium (Blöthe and Korup 2013). This geometry can be estimated based on borehole data that reports alluvium thickness in valley bottoms, employing geostatistical methods such as kriging (e.g., Deleplancque et al. 2018).

Kriging techniques are widely employed in geosciences to ensure the spatial continuity of properties beyond referenced observations (Matheron 1962, 1963; Chilès and Desassis 2018). These techniques interpolate the value of a variable at a given location by applying a weighted sum of observations of the variables at known locations. The weights account for the spatial structure of the data and the relative positions of the observations and the target point. Kriging has been successfully used to reconstruct alluvial bottom geometry from borehole data within the Seine River catchment, North France (Labarthe 2016; Deleplancque et al. 2018). One limitation of stationary kriging is the need to adjust correlation ranges for a fixed regional direction. While this approach is suitable for nearly straight valleys (Deleplancque et al. 2018), it becomes problematic in sinuous valleys where the direction continuously changes, particularly in elongated and narrow valleys. The use of regional variograms along a linear axis of anisotropy can lead to linear features (Seifert and Jensen 2000). Consequently, the predicted alluvial thickness may exhibit an unrealistic linear geometry due to discrepancies between the defined regional linear axis and the local direction of the river course (Labarthe 2016). Moreover, this direction may be unsuitable for estimating the kriging parameters such as the maximum and minimum ranges. Instead, an ideal kriging pattern should align with the average paleo-flow direction that shaped the valley.
In this study, we compare kriging methods for reconstructing alluvial bottom geometries using stationary and nonstationary anisotropy fields. The methods are tested along the Oise valley, a sinuous valley with frequent variations in bedrock lithology, located in the northern France. The resulting (i) bottom valley alluvium thickness and (ii) bedrock elevation maps are presented, along with their associated uncertainties, and discussed in relation to the Late Quaternary river dynamics of the study area.

## 2 Study Area

The study area encompasses the Oise valley, situated in northern France (Fig. 1). Originating south of Chimay in Belgium, the Oise River is a 341 km long tributary of the Seine River. Our investigation focuses on an approximately 180 km-long segment of the downstream part of the Oise River valley, from its confluence with the Noirrieu River near Vadencourt to its junction with the Seine River in Conflans-Saint-Honorine (Fig. 1a). The Oise River flows through the Paris Basin, an intracratonic sedimentary basin that underwent gradual subsidence from Triassic to Cretaceous periods (Guillocheau et al. 2000; Briais et al. 2016) ending with the extensive deposition of hundreds of meters of Upper Cretaceous chalk (Lasseur 2007) that locally outcrops in the study area (Fig. 1b). A very low subsidence phase began from lower Paleocene until late Eocene and the activation of NW-SE structures (e.g. Vigny, Bray and Margny anticlines, Artois bulge; Fig. 1a) interpreted as a response to the Pyrenean collision (Guillocheau et al. 2000; Briais et al. 2016). Several Paleogene episodes of marine transgression and regression resulted in the deposition of alternating layers of conglomerates, sand and clays as well as limestones of either continental or marine

<figure><img src="image_3.png" /><figcaption>**Fig. 1** a Geological settings of the study area from the 1/1,000,000 geological map of France (Chantraine et al. 2003). The locations of boreholes (n = 1010) are shown with black dots. **b** Simplified sedimentary log of the Oise valley bedrock lithologies (after Maton et al. 2014)</figcaption></figure>
origin within the Paris Basin (Pomerol 1989; Huyghe et al. 2015; Briais et al. 2016; Moreau et al. 2024). These alternations are observed in the Oise valley above the Creaceous chalk (Fig. 1). From the lower Oligocene to the Miocene, the subsidence is then reoriented along NE-SW axes while NE-SW flexure or buckling is interpreted as consequences of the Alpine collision (Guillocheau et al. 2000; Bourgeois et al. 2007).

Quaternary incision dynamics are largely influenced by climatic variations and moderate uplift in the study area (Larue 2000; Cojan et al. 2007; Antoine et al. 2007). A long-term incision rate of 50–60 m/Ma is inferred from chronology of valley slope terraces, with fluvial erosion occurring primarily at the transition from interglacial to glacial periods (Antoine et al. 2007). In the Oise Valley, existing radiocarbon and Optically Stimulated Luminescence (OSL) dating in the valley bottom is sparse, with the oldest known alluvium dating back to 32,750 years BP (Valladas 1981, 1994; Audouze et al. 1981; Pastre et al. 1997, 2000, 2002; Larue 2000; Granai 2014; Blaser et al. 2017; Granai and Limondin-Lozouet 2018; Table 1). Consequently, it is likely that valley bottom geometries and alluvium infill are largely inherited from the incision during the last glacial period and subsequent aggradation during the Late Glacial and Holocene periods (Pastre et al. 1997, 2000; Krier 2004; see Fig. 2). Recent OSL dating of fluvial sands in Maurecourt, located in the lower part of the Oise Valley, returned an age corresponding to Marine Isotope Stage 2 (MIS 2; Blaser et al. 2015; Chourio Camacho 2024), suggesting that Holocene deposits may be less widespread than originally thought by Krier (2004).

# 3 Methods

Kriging predictions and geostatistical simulations were programmed with Python scripts that called the gstlearn library (MINES Paris - PSL University 2025). The Jupyter Notebook files used in the study and the data are available on a website (see Acknowledgments).

## 3.1 Datasets

Kriging was performed using two types of data: (i) a polygon representing the limits of the alluvial plain (Fig. 3a), and (ii) a set of 1010 points containing information about the alluvium thickness (Fig. 3b).

The limits of the alluvial plains were extracted from vectorized and harmonized 1/50,000 scale geological maps of the French geological survey (BRGM). The youngest alluvial deposit polygons (i.e. *Fx*, *Fy* and *Fz* terraces) were merged to form the valley bottom polygon using the ArcGIS® software. Additionally, terraces which are more than 10 m above the river were not included within the valley bottom polygon to limit the study to the latest evolution stage of the valley. Finally, an inconsistency was corrected to the South of the Thieux anticline, near Verberie (Fig. 1) where the Oise valley is about 5 km wide. There, sands from the Ypresian (i.e., Bracheux sands) were mistakenly reported as alluvial deposits on the geological map (Compiègne and

<figure><img src="image_4.png" /></figure>
<table><caption>**Table 1** Synthesis of existing ages in the Oise alluvial plain</caption><thead><tr><th>Location</th><th>Dating methods</th><th>Dated material</th><th>Age (ka)</th><th>Error (ka)</th><th>References</th></tr></thead><tbody><tr><td>Boran-sur-Oise</td><td>radiocarbon</td><td>Wood stack</td><td>3.88</td><td>0.07</td><td>Krier 2004</td></tr><tr><td>l'Isle-Adam</td><td>OSL</td><td>Fine sand</td><td>10.28</td><td>0.80</td><td>Blaser et al. 2017</td></tr><tr><td>Croix-Saint-Ouen</td><td>radiocarbon</td><td>Organic sediment</td><td>12.40</td><td>0.12</td><td>Pastre et al. 2000</td></tr><tr><td>Maurecourt</td><td>OSL</td><td>Fine sand</td><td>17.00</td><td>1.2</td><td>Chourio Camacho 2024</td></tr><tr><td>Choisy-au-Bac</td><td>radiocarbon</td><td>Organic sediment</td><td>5.12</td><td>0.04</td><td>Granai and Limondin-Lozouet 2018; Granai 2014</td></tr><tr><td>Verberie</td><td>radiocarbon</td><td>Organic sediment</td><td>12.50</td><td>0.01</td><td>Pastre et al. 1997 (Audouze et al. 1981; Valladas 1981, 1994)</td></tr><tr><td></td><td>thermoluminescence</td><td>Silty sand</td><td>13.00</td><td>1</td><td></td></tr><tr><td>Béthisy-Saint-Martin</td><td>radiocarbon</td><td>Organic sediment</td><td>11.14</td><td>0.10</td><td>Larue 2000</td></tr><tr><td>Bruyères</td><td>radiocarbon</td><td>Organic sediment</td><td>4.25</td><td>0.28</td><td></td></tr><tr><td></td><td></td><td></td><td>10.28</td><td>0.11</td><td></td></tr><tr><td>Les Epinières Houdancourt</td><td>radiocarbon</td><td>Plant remains</td><td>12.54</td><td>0.11</td><td>Pastre et al. 2002; Ponel et al. 2005</td></tr><tr><td></td><td></td><td>Wood</td><td>12.06</td><td>0.11</td><td></td></tr><tr><td></td><td></td><td>Wood</td><td>11.62</td><td>0.09</td><td></td></tr><tr><td></td><td></td><td>Plant remains</td><td>11.26</td><td>0.09</td><td></td></tr><tr><td></td><td></td><td>Peat</td><td>10.90</td><td>0.14</td><td></td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Cross-sections in the alluvium of the lower part of the Oise valley (after Krier 2004). See Fig. 1 for location. Based on existing and recent dating (Blaser et al. 2015; Chourio Camacho 2024) of alluvial sediments (Table 1), the extension of Holocene sediment has been largely reduced compared to the initial interpretation by Krier (2004)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3** Datasets used to perform directional kriging predictions and simulations. **a** Data points, colored by elevation of the bedrock-alluvium interface, and polygon of the limits of the alluvial plain (see text for details). **b** Histogram of alluvium thickness picked from the boreholes. **c** Zoom in the Verberie area showing the construction of the anisotropy field based on the contours of the alluvial plain. Note that a 1000 m resolution is used for visualization while it is 100 m in the simulations</figcaption></figure>

<figure><img src="image_7.png" /></figure>
Senlis; http://infoterre.brgm.fr/) instead of valley bottom alluvium as shown by Pastre et al. (2000) and confirmed by our field observations and analyses of the boreholes.

Borehole data were extracted from the BSS (Banque du Sous-Sol) database available on the French Geological Survey website (http://infoterre.brgm.fr/). Each borehole was manually analyzed to ensure data consistency. Inconsistent data, accounting for approximately 20% of the analyzed boreholes, were discarded.

## 3.2 Construction of the Orientation Field

To perform the nonstationary kriging, an orientation field was generated based on the limits of the alluvial plain, based on the hypothesis that these limits are consistent with the paleo-river flow that carved the valley. This polygon was first used to build a centerline of the alluvial plain using the ArcGIS® software and was then split into two polylines (i.e. northern and southern edges of the plain; Fig. 3c). Two additional lines of control points bordering the outside of the plain to the North and to the South (with a 45° N direction) were also added. As a result, five polylines were constructed. Tangential vectors at each vertex of these polylines were generated using the NumPy® gradient method and then converted into angles using the NumPy® arctan2 function. Finally, these angles were interpolated linearly into a structured map using the SciPy method for interpolating unstructured data to grid (Fig. 3c).

## 3.3 Kriging Predictions

### 3.3.1 Principle of the Methods

Kriging is an interpolation method widely used in spatial statistics to predict the value of a variable defined across a spatial domain using a finite set of point observations (Wackernagel 2003). More precisely, the value of the variable at an unobserved location is predicted using a weighted sum of the observations, with weights computed to ensure that the resulting predictor is unbiased and has minimal variance. In particular, these weights account for both the (spatial) correlation of the variable between the observed locations, and the correlation between the prediction locations and the observed locations (Chilès and Desassis 2018).

In covariance-based kriging, one specifies a spatial correlation function, which gives the correlation of the variable between any two points of the spatial domain, thus allowing to compute the weights described above for any prediction location in the domain (Fig. 4). A usual assumption, called isotropy, is to consider that this function only depends on the distance between the points. If the correlation also depends on the orientation between the two points, the phenomenon is said to be anisotropic (Fig. 4a). In both cases, a variogram analysis can be used to straightforwardly estimate the correlation function from the observations (Wackernagel 2003).

However, in most real cases, stationarity is not a realistic assumption, as the spatial structure of the variable is expected to change across the domain. For instance, the correlation of the variable may be higher along a given direction, which might not be the same depending on where we are in the domain. This is the case in our work
<figure><img src="image_8.png" /><figcaption>Fig. 4 Correlation map of the 3 kriging methods tested in this study. a Anisotropy fields and b resulting correlation for the different models</figcaption></figure>

where we expect these directions to be defined by the vector field, which follows the valley, introduced in the previous subsection (Fig. 3). To account for these so-called local anisotropies, several approaches exist.

In this study, three spatial methods were investigated and compared (Figs. 4 and 5). The three methods used a Matérn covariance function with smoothness 1 for easy comparison with each other's. For the sake of the demonstration, kriging methods were performed with or without a drift order of 1, i.e., with latitude and longitude used as external covariates (Table 2). To estimate the intercept (mean) parameter of the drift, the coordinates have been centered to the middle of the domain. The borders of the domain were extended to reduce side effects. An uncorrelated (i.e., pure nugget) model was also used for comparison.

In the first approach (Cov-based stationary), the anisotropy of the prediction is considered stationary (Fig. 4a). In the second one (Cov-based nonstationary), the

<figure><img src="image_9.png" /></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 5** Kriging maps predicting alluvium thickness. **a** Covariance-based kriging with stationary geometric anisotropy. The estimated constant angle of anisotropy is 23.8°. **b** Covariance-based kriging with non-stationary geometric anisotropies. **c** SPDE kriging with non-stationary geometric anisotropies. **d** Zoom of the different kriging maps (**a**–**c**) in the Cergy-Pontoise area (meandering type). **e** Zoom of the different kriging maps (**a**–**c**) in the Verberie area (braided type)</figcaption></figure>

covariance is locally stationary, but the anisotropy can vary along the field. The covariance between the two points is then computed from the local models at these two points (Paciorek and Schervish 2006; Fouedjio et al. 2015). The last approach is based on Stochastic Partial Differential Equations (SPDE nonstationary) and parametrized similarly to the previous nonstationary one (Fig. 4b) but the covariance between two points results from the model over the field (i.e., along a path that follows the anisotropy field). This approach does not require to specify a global correlation function for the variable, nor the use of moving neighborhoods. Instead, the variable is modeled as (a sample of) the solution of a specific SPDE (Lindgren et al. 2011). The parameters of this SPDE, which can vary across the domain to account for local anisotropies, are in turn directly linked to the local correlation structure of the variable (Fuglstad et al. 2015). In practice, closed-form formulas allow to directly compute (inverse of) covariance matrices of the variable, which in turn are used to compute the same kriging weights as before (Pereira et al. 2022). In this study, we use triangulation and finite element method approximations for SPDE calculations. We did several tests by increasing the resolution of the mesh, until no significant improvement was observed. The resulting
<table><caption>**Table 2** Summary of parameters in the different kriging methods</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="3">Models without drift</th><th colspan="3">Models with drift</th></tr><tr><th>Cov-based stat</th><th>Cov-based nonstat</th><th>SPDE nonstat</th><th>Uncorrelated</th><th>Cov-based stat</th><th>Cov-based nonstat</th><th>SPDE nonstat</th><th>Uncorrelated</th></tr></thead><tbody><tr><td rowspan="10">Estimated drift parameters</td><td>Mean (m)</td><td>6.91</td><td>6.89</td><td>6.90</td><td>7.24</td><td>6.73</td><td>6.70</td><td>6.71</td><td>6.52</td></tr><tr><td>Estimated covariance function parameters</td><td>X</td><td>x</td><td>x</td><td>x</td><td>x</td><td>5.42E-07</td><td>4.28E-06</td><td>2.72E-06</td><td>5.10E-05</td></tr><tr><td></td><td>Y</td><td>x</td><td>x</td><td>x</td><td>x</td><td>-2.85E-05</td><td>-3.32E-05</td><td>-3.10E-05</td><td>-1.03E-04</td></tr><tr><td>Angle (°)</td><td>63.4</td><td>x</td><td>x</td><td>x</td><td>x</td><td>64.7</td><td>x</td><td>x</td><td>x</td></tr><tr><td>Range 1 (m)</td><td>1470</td><td>1972</td><td>1688</td><td>x</td><td>x</td><td>1230</td><td>1631</td><td>1436</td><td>x</td></tr><tr><td>Range 2 (m)</td><td>985</td><td>821</td><td>858</td><td>x</td><td>x</td><td>887</td><td>735</td><td>763</td><td>x</td></tr><tr><td>Range ratio</td><td>1.49</td><td>2.40</td><td>1.97</td><td>x</td><td>x</td><td>1.39</td><td>2.22</td><td>1.88</td><td>x</td></tr><tr><td>Variance</td><td>5.34</td><td>5.35</td><td>5.43</td><td>x</td><td>x</td><td>4.83</td><td>4.86</td><td>4.94</td><td>x</td></tr><tr><td>Variance error</td><td>2.81</td><td>2.82</td><td>2.84</td><td>8.34</td><td>x</td><td>2.77</td><td>2.77</td><td>2.79</td><td>7.46</td></tr><tr><td>Total Variance</td><td>8.15</td><td>8.17</td><td>8.27</td><td>8.34</td><td>x</td><td>7.60</td><td>7.63</td><td>7.73</td><td>7.46</td></tr><tr><td rowspan="4">Likelihood</td><td>MAE</td><td>-2334.1</td><td>-2329.4</td><td>-2330.5</td><td>-2541.4</td><td>-2325.6</td><td>-2321.3</td><td>-2322.5</td><td>-2484.3</td></tr><tr><td>Cross-validation (random)</td><td>1.67</td><td>1.65</td><td>1.65</td><td>2.30</td><td>1.68</td><td>1.65</td><td>1.66</td><td>2.20</td></tr><tr><td>RMSE</td><td>2.23</td><td>2.22</td><td>2.22</td><td>2.89</td><td>2.24</td><td>2.23</td><td>2.23</td><td>2.74</td></tr><tr><td>Cross-validation (spatial)</td><td>2.09</td><td>2.05</td><td>2.06</td><td>2.31</td><td>2.05</td><td>2.02</td><td>2.02</td><td>2.22</td></tr><tr><td></td><td>RMSE</td><td>2.73</td><td>2.69</td><td>2.69</td><td>2.90</td><td>2.63</td><td>2.61</td><td>2.61</td><td>2.76</td></tr></tbody></table>

<figure><img src="image_11.png" /></figure>
mesh has 203,428 vertices. The resolution is 100 m along the valley and 50 m in the
orthogonal direction.

3.3.2 Estimation of the Ranges and Kriging Maps

To perform the kriging (Fig. 5), a matrix is necessary to describe the anisotropy. It corresponds to a 2D ellipsoid where the semi-major and semi-minor axes are called the maximum and minimum ranges. These ranges were estimated before the kriging by (profiled) maximum likelihood under gaussian assumption (Diggle and Ribeiro 2007). We decided not to include the borders of the valley to infer these parameters as they would induce bias in both the parameter estimation and the calculation of volumes. For the covariance-based stationary case (Fig. 5a), the main orientation is an additional parameter estimated by maximum profiled likelihood. For the nonstationary cases (covariance-based and SPDE) (Figs. 5b and 5c), the predictions were made using the orientation (angles) field (Table 2). For the stationary case, the anisotropy angle was constant.

3.3.3 Cross-Validation and Assessment of Uncertainties

To quantitatively compare the outputs of the different methods and to assess their
predictability, k-fold cross-validation tests were performed. Each cross-validation test
was repeated around a hundred times for each method. For each kriging method, a
subset was selected in two ways: (a) by removing randomly 1% of the data or (b) by
removing data spatially inside a sliding window corresponding to 1% of the length
of the valley. This latter method allowed testing the predictability of kriging without
neighboring constraints. In the two cases, parameters were inferred on the whole
dataset once for all. Prediction ability of the methods at the location of the removed
data points was assessed by calculating either the root mean square error (RMSE) or
the mean absolute error (MAE; Table 2).

3.3.4 Alluvium-Bedrock Interface

The bedrock-alluvium interface elevation was then reconstructed by subtracting the kriging prediction of alluvial thickness to the elevation of the top of the alluvial plain extracted from a 500 m resolution DEM (i.e., degraded STRM 3 arcsecond resolution) using the Raster to Point routine in ArcGIS®. A SPDE kriging was performed on the 2470 evenly spaced points extracted from the DEM. The Topographic Analysis Toolbox (TAK; Forte and Whipple 2019, based on the TopoToolbox Matlab library of Schwanghart and Kuhn 2010; Schwanghart and Scherler 2014) was used to extract swath profiles along the Oise Valley (Fig. 6). Only the minimum elevations were kept as the thalwegs.

3.3.5 Alluvium Volume

The alluvial volume was estimated based on a set of 500 SPDE nonstationary simu-
lations applying the inferred parameters (Fig. 7).
<figure><img src="image_12.png" /></figure>

**Fig. 6** Elevation (SWATH) profiles along the centerline of the Oise valley showing the top of the alluvial plain (in grey) and the valley thalwegs at the bedrock-alluvium interface computed with the TAK (Forte and Whipple 2019). Green rectangles correspond to areas with chalk bedrock lithology. The major confluences are indicated by blue arrows

<figure><img src="image_13.png" /></figure>

**Fig. 7** Justification of the approach (a, b) and calculation (c) of the volumes performed based on 500 SPDE conditional simulations of the same model. a, b Sketches of a section with boreholes illustrating the difference between kriging with or without borders. The mean (i.e., the value towards which the kriging will tend when no other data constraint is available) estimated with borders (where thickness equals zero) is biased to a lower value compared to the mean estimated without border. This does not affect too much the kriging when data is available, although some edge effects may occur (a). In the section with only borders, the result is biased twice because only zero values are modelled and because the mean is lower. Kriging was therefore performed without borders in this study

<figure><img src="image_14.png" /></figure>
# 4 Results and Interpretations

## 4.1 Parametric Analysis

### 4.1.1 Covariance Function Parameters

The results of the covariance matrix estimations are shown in Table 2. For the models with drift, the log likelihood parameter is higher for the covariance-based and SPDE (−2321.3 and −2322.6) nonstationary krigings compared to the stationary one (−2325.6). According to these values, the best parameters are estimated using non-stationary predictions. Variance, variance error and total variance are comparable in all cases, although a bit higher for SPDE (4.94, 2.79 and 7.73 respectively). Krigings without drift have a systemically lower likelihood than their equivalent models with drift (Table 2).

Overall, the models can be ranked by order of pertinence based on their log likelihood parameter values. First, models with spatial correlation are better than uncorrelated models. Second, models with drifts are systematically better than models without drift. Finally, non-stationary models are better than the stationary one.

The predicted ranges are between 1230 and 1631 m for the maximum one and between 735 and 887 m for the minimum one (Table 2). The ranges ratio is smaller in the stationary case (i.e.,1.4 along the 23.8° direction) with a maximum range of 1230 m and a minimum range of 887 m. By comparison, a stationary estimation of these ranges based on a variogram fitted with a Matérn covariance function returned values of 896 m and 1204 m along the 24° and 114° directions, respectively. With the non-stationary methods the estimated maximum ranges were 1631 m and 1435 m for the Cov-based and SPDE methods, respectively returning ranges ratios of 2.22 and 1.88. In covariance-based approach the anisotropy ratio thus is higher. It is explained by the fact that when there is a change of anisotropy difference between two points, continuity is “broken” between the two points in covariance-based approach while there is more continuity in SPDE (Fig. 4). Finally, krigings without drift have a higher range ratio (Table 2). This can be explained by the fact that the drift captures some of the anisotropy, which is compensated by increasing ranges and range ratio in models without drift.

### 4.1.2 Cross-Validation

Results of the cross-validations (MAE and RMSE) show that the kriging methods using non-stationary anisotropies are better predictors, both randomly and spatially (Table 2). The MAE and RMSE values are indeed higher with stationary anisotropies while they are lower and similar for the non-stationary cases (covariance-based and SPDE). Finally, and as expected, the prediction without spatial correlation has the highest errors.

Cross-validations indicate that the use of the drift improves the results in the spatial case and not in the random case. Indeed, in the latter case, observations are available everywhere and mean specifications have only a small influence whereas the spatial case tends to mimic extrapolation, where drift compensates for the lack of information.
## 4.2 Geometries of the Kriging Maps

In all predictions, the linear drift accounts for the increase of alluvial thickness from upstream to downstream (Fig. 5). The map of alluvium thickness generated with the covariance-based kriging with stationary anisotropies shows an overall patchy geometry (Fig. 5a). A slight anisotropy is observed along the 24° direction, consistent with the calculated ranges. In some cases, the predicted geometries are very rounded, patchy, such as in the Verberie, Compiègne-Noyon or Chauny areas (Figs. 1 and 5a). In cases where the density of data is low, the patchiness is increased. Overall, the pattern – although consistent with the data distribution – lacks continuity and produces geometries that are not expected for a river system flowing into a valley and carving into the bedrock.

The nonstationary predictions exhibit more continuity along features that follow the valley direction, compared to the kriging with stationary anisotropies. It is very visible in the Verberie area, where the kriging with stationary anisotropies has a very regular pattern, whereas the two kriging maps with non-stationary anisotropies mimic the continuity of the river path that carved the valley (Fig. 5d). In the curved valley of the Cergy-Pontoise area, continuity is also improved in the two non-stationary kriging maps, although more subtly (Fig. 5e). Slight differences exist between the covariance based and SPDE predictions with non-stationary anisotropies. Features are indeed narrower and more elongated in the covariance-based kriging. It is due to the differences in estimated ranges.

As shown above, the pattern of alluvial thickness of the classical kriging is rather uniform and patchy. In krigings with nonstationary anisotropies, distinct types of geometry are observed. The first type is visible in wide alluvial plains, particularly near Verberie (Figs. 5d and 5e). There, alluvial thickness varies across the width of the valley. In the short-range cases, it mimics a braided pattern (Fig. 5e). The width and length of the associated bars are in the order of 750–1500 m and 1000–3000 m respectively. Similar patterns are observed near Beaumont-sur-Oise and upstream of Chauny (Fig. 5b and c). The other type is the alternation of topographic highs and lows along narrow, sinuous valley sections that have a riffle and pool geometry (Wolman and Leopold 1957). This pattern is dominant downstream of Beaumont-sur-Oise, particularly in the Cergy-Pontoise area (Fig. 5d). There, the spacing between pools (i.e., topographic lows or thickness highs) is 3000–7000 m. Similar patterns are observed near Compiègne and Creil (Fig. 5b and c).

The Oise Valley thalwegs (Fig. 6) exhibit contrasted geometries, consistent with the observations made on the maps. Different segments of the valley can be observed. The upstream section of the Oise Valley (i.e., > 150 km) has a steeper slope and smaller alluvial thickness compared to the downstream section (i.e., 0–150 km) (Fig. 7). This corresponds to the transition to chalk bedrock but also to the confluence with the Serre River. Finally, the two identified styles (braided and meandering) can be distinguished from the thalweg: i) riffle-pool patterns are visible with ups and downs (up to 10 m in elevation) of the thalweg, mostly from Cergy to Beaumont-sur-Oise but also near
Compiègne, while ii) braided-like patterns correspond to sections with thalweg elevation oscillations of shorter amplitude (less than 5 m), visible in the Verberie and Chauny areas.

Finally, the volume estimated by simulations is about $2.30 \pm 0.02 \text{ km}^3$ (Fig. 7). By comparison, the volume calculated from the kriging map is $2.34 \text{ km}^3$. As explained in Fig. 7 and developed in the discussion, it is considered as a maximum estimation.

# 5 Discussion

## 5.1 Relevance of the Kriging with Non-stationary Anisotropies to Predict Alluvial Thickness

In this study, kriging using nonstationary geometric anisotropies appears adapted to predict alluvial thicknesses geometries in the case of the moderately sinuous Oise valley. Moreover, the vector field used, which is based on geological contours of the valley bottom, is relevant to improving the kriging. As shown in Table 2, non-stationary hypotheses allowed inferring sets of parameters with higher likelihood. The inferred ranges ratios are higher in nonstationary cases than in the stationary one. This supports that the anisotropy field is relevant to better capturing the anisotropies and to improving the kriging. Also, the cross-validation indicates that the two non-stationary methods (SPDE and covariance-based) have systematically less uncertainties than the stationary one, although the difference is small (2–5 cm). One should note that there are no major differences between cross-validations of the two nonstationary methods (Table 2).

Estimated anisotropy range ratios are lower in the stationary case (1.4) compared to the non-stationary cases (1.88 and 2.22 for SPDE and Cov-based cases, respectively). These anisotropy ranges are consistent with known values in natural systems, although a bit low. Indeed, for braided rivers, the average aspect ratio of bars (length divided by width) is about 5 (Kelly 2006), although some dispersion is observed, with values from 2 to 10 (Limaye et al. 2018), while the confluence-bifurcation spacing is 4–5 times the width of channels (Hundey and Ashmore 2009). For meandering rivers, recent compilations indicate that the ratio of the meander wavelength to the bankfull width of river cross-section is ca. 11 (Lemay et al. 2020). In the Beaumont-sur-Oise area, the high density of boreholes has allowed for the construction of very detailed cross-sections of the alluvial bottom, particularly the geometry of the bedrock-interface (Krier 2004). Locally, very densely spaced boreholes (along cross-sections realized in the frame of subway tunnels) indicate rapid fluctuations of the elevation of the erosional bedrock-alluvium interface. There, the width and spacing of channelized erosion features that can be individualized is in the order of 200–1000 m (Fig. 2), indicating the relevance of the minimum ranges found (700–900 m).

Due essentially to their high range ratio and the orientation of the field of angles, the kriging maps with nonstationary anisotropies have a wealth of detailed and realistic geometries with regards to river geomorphology knowledge (Fig. 5). For instance, one can observe braided-like patterns along wide valleys, and riffle-pool patterns along narrow and sinuous valley bottoms, consistently with the existing theory (e.g. Parker
1976). In the former, the resulting aspect ratio of “braids” (3–5) fits well with observations in natural systems (Kelly 2006). Finally, alluvial thickness increases immediately after the confluences with the Ailette and Serre rivers, likely reflecting scouring, as observed in the field (Best and Ashworth 1997) and experimentally (Grimaud et al. 2017). Such deepening is also observed with the geometry of the bedrock at the confluence with the Aisne River (Fig. 7). Together, these observations point towards the specifics of flow and sediment transport carving the bedrock at the valley floor.

Hence, the kriging with nonstationary spatial anisotropies introduced in this study has shown its relevance to predicting alluvium thickness along sinuous valley bottoms. Although the reduction of uncertainties (shown with the cross-validation) is modest, improved spatial realism is potentially very relevant for hydraulic simulations, particularly because the model anisotropies may affect flow path in the alluvium. Another interesting point is that patterns consistent with fluvial geomorphology, with local differences (i.e., there are meandering vs braided-dominated areas) emerged directly from a regional set of parameters. Much remains to be done in the future to improve this approach. As shown by the literature, ranges and range ratio may be different when the fluvial style changes. Furthermore, there is a tendency for bars and scours size to increase with drainage area. Hence, instead of a constant set of maximum and minimum ranges, it could be possible to infer sets of ranges that evolve (increase) with drainage area. Finally, the anisotropy field used in this study was obtained from geological mapping. Such data, although reasonably well mapped, may be locally inaccurate. In the future, it may be interesting to compute the anisotropy field directly from the data, and to compare it with one based on geological mapping.

## 5.2 Geological Implications

The geometries of alluvial deposits observed in the Oise valley bottom likely reflect paleo-hydraulic conditions from Glacial to Late Glacial times, the Meghalayen representing only limited modifications of the alluvial plain morphologies in the Seine River catchment (Pastre et al. 1997; 2000; Antoine 1997). This is supported by the dating of the deeper overlying sands, which belong to MIS2 period in the downstream Oise valley in Maurecourt (Fig. 3). The braided patterns observed in the non-stationary kriging maps are consistent with the interpretation that fluvial systems of Western Europe, the Seine River system in particular, were dominated by braided-like systems (Bridgland and Westaway 2008; Antoine et al. 2007). These morphologies observed at the alluvium-bedrock interface would correspond to erosion inherited from glacial periods. A novelty is the observed riffle-pool geometry at the scale of the valley along narrow valley sections. Such geometry is well known for rivers (Richards 1976) but not so much at the scale of valleys. It is illustrated by the geometry of the valley thalweg, which displays ups and downs that are saliant in the downstream part of the Oise valley near Cergy and Pontoise (Fig. 7). This supports a link between valley aspect ratio and valley bedrock patterns (i.e., braided vs meandering), as shown for rivers (Parker 1976). Ultimately, it points out that braided patterns may not be systematically associated with glacial deposits, an assertion that must be verified with the facies analysis of alluvial deposits, in the valley bottom or alluvial terraces.
At the scale of the valley, the change in slope and alluvial thickness upstream of Chauny (Fig. 7) confirms the coincidence of bedrock lithology change and alluvial/bedrock slope change (Grimaud et al. 2014). This transition also corresponds to the confluence with the Serre River, and there could be an equally important effect of increasing discharge after the confluence on the observed changes (i.e., Crosby et al. 2007). Krier (2004) speculated that the thalweg in the downstream Oise valley would be “over”-incised and would mark regressive incision attributed to condition before the Late Glacial. In the Cergy-Pontoise, the thalweg profile—perhaps limited by the amount of available borehole data—indicates that erosion did not erase the riffle-pool geometry. Unfortunately, the riffle-pool geometry makes it difficult to accurately identify knickpoints that would be related to regressive erosion in response to base level drop during glacial periods, as suggested by Deleplancque et al. (2018) along the Bassée alluvial plain (Fig. 7). Alluvial thickness rarely exceeds 10 m upstream of Compiègne and the few areas with increased thickness are related to confluences (i.e., Serre and Ailette Rivers). One may thus tentatively suggest that a transition exists there, but it remains speculative. Narrowing and shallowing of the valley bottom was shown upstream of Beaumont-sur-Oise by Larue (2000) who suggested that this could be the result of Late Quaternary uplift of the Bray anticline. The kriging maps show geometry that is consistent with this pattern. However, it also shows other ups and downs of the bedrock elevation elsewhere, which are related to meandering pattern. Although enhanced when compared to other areas, the increased elevation of the bedrock when crossing the Bray anticline is not inconsistent with the active uplift but may also simply be due to the local change of bedrock lithology to chalk as opposed to Cenozoic sediments (Fig. 7).

Finally, the kriging maps and estimations allow calculating the volume of alluvium temporarily stored in the Oise valley. As shown by Fig. 6, the output is a maximum compared to one computed with the borders of the plain used as zero values for alluvium thickness. However, including the zeros leads to many biases, especially in sections with only the borders as constraints. There, only zero values are modelled, and the mean -towards which the kriging will tend without further constraints from data- is lower. This results in an underestimation of the volume. Kriging was therefore performed without borders. The estimated 2.3 km³ volume of alluvium is therefore the best proxy. It has a similar bias towards the mean when no data is available, but this bias is weaker. Knowing this volume is of importance to estimate the temporary storage in transfer zones from a source-to-sink point of view (Castelltort and van den Driessche, 2003; Blöthe and Korup 2013), but also to constructing, gravel-mining, water supply companies, and for public authorities who manage and supervise valley bottom landscapes.

# 6 Conclusion

In this study, we test different kriging methods for predicting alluvial thickness geometries in the valley bottom using random fields with nonstationary geometric anisotropies constructed from the contours of the alluvial plain. The two methods based on covariance and Stochastic Partial Differential Equations are compared to
kriging with stationary anisotropy. The main parameters are estimated by maximum likelihood and the different methods are tested using cross-validations. The methods with nonstationary anisotropies present fewer uncertainties than the one with stationary anisotropies and produces patterns realistic with river geomorphic and sedimentologic knowledge. Evidence of braided patterns, riffle-pool geometry related to meandering, and scouring at confluences is observed at the scale of the valley width. Finally, it is estimated that about 2.30 km³ of alluvium is currently stored in the studied Oise valley bottom. Although the method is applied to an alluvial example, it could easily be adapted to other geological contexts where anisotropy field can be estimated, such as faulted and folded structures, supergene mineralization or environmental context to track pollutants dissemination in rivers.

**Acknowledgements** This work is part of the RGF: Bassin de Paris program supported by the BRGM. It benefited from financial support from the Agence de l'Eau Seine Normandie. We thank Alexandre Leveque and Isabelle Cojan for fruitful exchanges early in this work. gstlearn version 1.7.0 was used in this study (https://github.com/gstlearn/gstlearn). A copy of the notebook used in this study can be found here: https://soft.mines-paristech.fr/gstlearn/oise_SPDE/

**Funding** Open access funding provided by Mines Paris - PSL.

**Declarations**

**Conflict of interest** All authors declare that they have no conflicts of interest.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

**References**

Antoine P (1997) Modifications des systèmes fluviatiles à la transition Pléniglaciaire-Tardiglaciaire et à l'Holocène: l'exemple du bassin de la Somme (Nord de la France). Géographie Phys Quat 51:93-106

Antoine P, Lozouet NL, Chaussé C, Lautridou J-P, Pastre J-F, Auguste P, Bahain J-J, Falguères C, Galehb B (2007) Pleistocene fluvial terraces from northern France (Seine, Yonne, Somme): synthesis, and new results from interglacial deposits. Quat Sci Rev 26:2701-2723

Audouze F, Cahen D, Keeley LH, Schmider B (1981) Le site magdalénien du buisson campin à verberie (Oise). Gall Préhistoire 24:99-143

Best JL, Ashworth PJ (1997) Scour in large braided rivers and the recognition of sequence stratigraphic boundaries. Nature 387:275-277

Best J, Darby SE (2020) The pace of human-induced change in large rivers: stresses, resilience, and vulnerability to extreme events. One Earth 2(6):510-514

Blaser F, Chaussée C, Raymond P, Blaser R (2015) Maurecourt «78 sentier des carreaux» (diagnostic INRAP No. D108832)

Blaser R, Chaussé C, Roncin O, Rochart X, Bayle G, Lefèvre A, Raymond P, Roque C (2017) L'Isle-Adam intersection des avenues Paul Thoureau & Chemin Vert (diagnostic INRAP No. D101013)

Blöthe JH, Korup O (2013) Millennial lag times in the Himalayan sediment routing system. Earth Planet Sci Lett 382:38-46
Bourgeois O, Ford M, Diraison M, de Veslud CLC, Gerbault M, Pik R, Ruby N, Bonnet S (2007) Separation of rifting and lithospheric folding signatures in the NW-Alpine foreland. Int J Earth Sci 96:1003–1031.
https://doi.org/10.1007/s00531-007-0202-2

Briais J, Guillocheau F, Lasseur E, Robin C, Châteauneuf J-J, Serrano O (2016) Response of a low-subsiding intracratonic basin to long wavelength deformations: the palaeocene–early eocene period in the Paris Basin. Solid Earth 7:205–228

Bridgland D, Westaway R (2008) Climatically controlled river terrace staircases: a worldwide quaternary phenomenon. Geomorphology 98:285–315

Castelltort S, Van Den Driessche J (2003) How plausible are high-frequency sediment supply-driven cycles in the stratigraphic record? Sediment Geol 157:3–13

Chantraine J, Autran A, Cavelier C (2003) Carte géologique de la France à 1/1 000 000, 6e édition révisée. Orléans BRGM

Chilès J-P, Desassis N (2018) Fifty years of kriging. Handbook of Mathematical Geosciences: Fifty Years of IAMG, pp 589–612

Chourio Camacho DN (2024) Late quaternary incision dynamics and valley bottom geometry within the Seine River catchment. PhD thesis Mines Paris, 306p

Clubb FJ, Mudd SM, Milodowski DT, Valters DA, Slater LJ, Hurst MD, Limaye AB (2017) Geomorphometric delineation of floodplains and terraces from objectively defined topographic thresholds. Earth Surf Dyn 5:369–385

Clubb FJ, Weir EF, Mudd SM (2022) Continuous measurements of valley floor width in mountainous landscapes. Earth Surf Dyn 10:437–456

Cojan I, Brulhet J, Corbonnois J, Devos A, Gargani J, Harmand D, Jaillet S, Jouve A, Laurain M, Lejeune O (2007) Morphologic evolution of eastern Paris Basin: «ancient surfaces» and quaternary incisions. Mém Société Géologique Fr 178:135–155

Crosby BT, Whipple KX, Gasparini NM, Wobus CW (2007) Formation of fluvial hanging valleys: theory and simulation. J Geophys Res Earth Surf 112:F03S10

Deleplancque B, Cojan I, Beucher H, Mehl C, Stab O (2018) Spatial and temporal patterns of the upper pleistocene alluvial fill deposits of the upstream Seine River alluvial plain, la Bassée, France. Geomorphology 318:148–161

Diggle PJ, Ribeiro PJ Jr (2007) Model-based geostatistics. In: Springer series in statistics. Springer New York

Forte AM, Whipple KX (2019) The topographic analysis kit (TAK) for TopoToolbox. Earth Surf Dyn 7:87–95

Fouedjio F, Desassis N, Romary T (2015) Estimation of space deformation model for non-stationary random functions. Spatial Statistics 13:45–61

Fuglstad G-A, Lindgren F, Simpson D, Rue H (2015) Exploring a new class of non-stationary spatial Gaussian random fields with varying local anisotropy. Stat. Sin. 115–133

Granai S, Limondin-Lozouet N (2018) The Holocene expansion of grassland in northern Europe reconstructed from molluscan assemblages. Boreas 47:768–779

Granai S (2014) L'anthropisation des milieux du Néolithique à l'âge du Fer dans le bassin de la Seine enregistrée par les malacofaunes continentales. UNIVERSITÉ PARIS 1 – PANTHÉON-SORBONNE

Grimaud J-L, Chardon D, Beauvais A (2014) Very long-term incision dynamics of big rivers. Earth Planet Sci Lett 405:74–84

Grimaud J, Paola C, Ellis C (2017) Competition between uplift and transverse sedimentation in an experimental delta. J Geophys Res Earth Surf 122:1339–1354

Grimaud J-L, Gouge P, Huyghe D, Petit C, Lestel L, Eschbach D, Lemay M, Catry J, Quaisse I, Emperor A (2023) Lateral river erosion impacts the preservation of Neolithic enclosures in alluvial plains. Sci Rep 13:16566

Guillocheau F, Robin C, Allemand P, Bourquin S, Brault N, Dromart G, Friedenberg R, Garcia J-P, Gaulier J-M, Gaumet F (2000) Meso-Cenozoic geodynamic evolution of the Paris Basin: 3D stratigraphic constraints. Geodin Acta 13:189–245

Hundey E, Ashmore P (2009) Length scale of braided river morphology. Water Resour Res 45(8):W08409

Huyghe D, Lartaud F, Emmanuel L, Merle D, Renard M (2015) Palaeogene climate evolution in the Paris Basin from oxygen stable isotope (δ18O) compositions of marine molluscs. J Geol Soc 172:576–587

Kelly S (2006) Scaling and hierarchy in braided rivers and their deposits: examples and implications for reservoir modelling. Blackwell Oxford, Oxford
Krier V (2004) La plaine alluviale de l'Oise: milieu et système fluvial; stratigraphie des formations alluviales et archéologie. Cons. Général Val D'Oise INRAP 2

Labarthe B (2016) Quantification des échanges nappe-rivière au sein de l'hydrosystème seine par modélisation multi-échelle

Langston AL, Temme AJ (2019) Impacts of lithologically controlled mechanisms on downstream bedrock valley widening. Geophys Res Lett 46:12056–12064

Larue J-P (2000) Morphodynamique fluviale et néotectonique dans la vallée de l'Oise (Bassin parisien, France). Bull Société Géologique Fr 171:577–585

Lasseur É (2007) La craie du Bassin de Paris (cénomanien-campanien, crétacé supérieur): sédimentologie de faciès, stratigraphie séquentielle et géométrie 3D. Rennes 1

Lemay M, Grimaud J-L, Cojan I, Rivoirard J, Ors F (2020) Geomorphic variability of submarine channelized systems along continental margins: Comparison with fluvial meandering channels. Mar Pet Geol 115:104295

Limaye AB, Grimaud J, Lai SY, Foreman BZ, Komatsu Y, Paola C (2018) Geometry and dynamics of braided channels and bars under experimental density currents. Sedimentology 65:1947–1972

Lindgren F, Rue H, Lindström J (2011) An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. J R Stat Soc Ser B Stat Methodol 73:423–498

Matheron G (1962) Traité de géostatistique appliquée. Editions Technip

Matheron G (1963) Traité de géostatistique appliquée. 2. le krigeage. Editions Technip

Maton D, Pannet P, Sandifort C (2014). Inventaire des mouvements de terrain du département de l'Oise et la Somme. (No. BRGM/RP-63956-FR), Rapport BRGM

MINES Paris - PSL University (2025) gstlearn/gstlearn: Stable Release v1.6.0 (stable_1.6.0). Zenodo. https://doi.org/10.5281/zenodo.14628813

Moreau K, Andrieu S, Briais J, Brigaud B, Ader M (2024) Facies distribution and depositional cycles in lacustrine and palustrine carbonates: the Lutetian-Aquitanian record in the Paris Basin. Depositional Rec 10:124–158. https://doi.org/10.1002/dep2.264

Paciorek CJ, Schervish MJ (2006) Spatial modelling using a new class of nonstationary covariance functions. Environ. off J Int Environ Soc 17(5):483–506

Parker G (1976) On the cause and characteristic scales of meandering and braiding in rivers. J Fluid Mech 76:457–480

Pastre J-F, Fontugne M, Kuzucuoğlu C, Leroyer C, Limondin-Lozouet N, Talon M, Tisnérat N (1997) L'évolution tardi et postglaciaire des lits fluviaux au nord-est de Paris (France). Relations avec les données paléoenvironnementales et l'impact anthropique sur les versants. Géomorphologie Relief Process Environ 3:291–312

Pastre J-F, Leroyer C, Limondin-Lozouet N, Chaussé C, Fontugne M, Gebhardt A, Hatté C, Krier V (2000) Le tardiglaciaire des fonds de vallée du bassin parisien (France). Quaternaire 11:107–122

Pastre J-F, Leroyer C, Limondin-Lozouet N, Orth P, Chaussé C, Fontugne M, Gauthier A, Kunesh S, Lejeune Y, Saad MC (2002) Variations paléoenvironnementales et paléohydrologie durant les 15 derniers millénaires: les réponses morphosédimentaire des vallées du Bassin Parisien (France)

Pereira M, Desassis N, Allard D (2022) Geostatistics for large datasets on riemannian manifolds: a matrix-free approach. J Data Sci 20(4):512–532

Pomerol C (1989) L'évolution du bassin parisien. assoc. sédimentologistes fr.-dyn. méthodes détude bassins sédimentaires Ed. Tech. Paris 165–178.

Ponel P, Coope R, Antoine P, Limondin-Lozouet N, Leroyer C, Munaut AV, Pastre JF, Guiter F (2005) Lateglacial palaeoenvironments and palaeoclimates from Conty and Houdancourt, northern France, reconstructed from Beetle remains. Quat Sci Rev 24(23–24):2449–2465

Richards K (1976) The morphology of riffle-pool sequences. Earth Surf Process 1:71–88

Schwanghart W, Kuhn NJ (2010) Topo Toolbox: a set of matlab functions for topographic analysis. Environ Model Softw 25:770–781

Schwanghart W, Scherler D (2014) TopoToolbox 2–MATLAB-based software for topographic analysis and modeling in earth surface sciences. Earth Surf Dyn 2:1–7

Seifert D, Jensen J (2000) Object and pixel-based reservoir modeling of a braided fluvial reservoir. Math Geol 32:581–603

Tofelde S, Bufe A, Turowski JM (2022) Hillslope sediment supply limits alluvial valley width. AGU Adv 3:e2021AV000641
Valladas H (1981) Thermoluminescence de grès de foyers préhistoriques: estimation de leur température de chauffe. *Archaeometry* 23:221-229

Valladas H (1994) Chronologie des sites du magdalénien final du bassin parisien: Le milieu naturel et son exploitation. I: climats et paysages. Doc. Archéologie Fr. 65-68

Wackernagel H (2003) Multivariate geostatistics: an introduction with applications. Springer Science & Business Media, New York

Wolman MG, Leopold LB (1957) River flood plains: some observations on their formation, in: US Government Printing Office Washington, DC. US Government Printing Office Washington, DC, pp. 87-109