Assessing lithological uncertainty in dikes: Simulating construction history and its implications for flood safety assessment

Dikes often have a long history of reinforcement, with each reinforcement adding new material resulting in a heterogeneous dike. As data on the dike internal heterogeneity is sparse, it is generally overlooked in the stability assessment of dikes. We present an object‐based and process‐based model simulating dike construction history on archeological dike cross, yielding similar patterns of heterogeneity as observed in real dikes, and apply it in a dike safety assessment. Model predictions improve when being based on more accurate statistics of dike buildup, or when being conditioned to ground truth data. When incorporated in a dike stability assessment, multiple model runs can be coupled to hydrological simulations and dike slope stability calculations, resulting in a probabilistic stability assessment considering internal dike heterogeneity. While high‐resolution observations are still sparse, good model accuracies can be reached by combining regional information on dike buildup with local point observations and this model provides a parsimonious basis to include information of internal dike heterogeneity in safety assessments.


| INTRODUCTION
Current and future changes in climate likely lead to more extreme precipitation events and high-water levels, which demand major investments in dike reinforcement projects. To efficiently guide these expenditures, extensive research is carried out to determine material properties of dikes (de Waal, 2016). This information is of great importance for dike stability, as the material type directly influences dike stability by the stability parameters (Cao et al., 2016), but also indirectly influences groundwater pressures and flow within the dike (Butera et al., 2020;Kuriqi et al., 2016;van Woerkom et al., 2021).
The material of which dikes are constructed is often very heterogeneous, as throughout their existence dikes were heightened and widened. To support this, a short overview of the river dike history of the Netherlands is presented, but a similar history can be expected elsewhere. The earliest dikes in the Netherlands were constructed in the early Middle Ages, and at the start of the 16th century all rivers were bordered by dikes. These dikes were owned by locals who could use the dike but were also responsible for its maintenance. A water board coordinated and inspected the maintenance works done by the local land owners (TeBrake, 2002). Large scale dike reinforcements only started in the 19th century with the mechanization in transport and construction and the establishment of a national institution for water safety (in Dutch: Bureau van de Waterstaat). In more recent times (following the 1953 storm surge) the Delta Law of 1958 led to an extensive round of dike reinforcement (Pleijster, Veeken and LOLA Landscape Architects, 2014). Currently, new risk-based safety standards based on acceptable probabilities of flooding for various dike failure mechanisms are defined, again leading to a series of reinforcement projects. Differences in reinforcement methods, and the prolonged lack of institutionalized control, has resulted in very heterogeneous dikes (Figure 1), and locally the original lithological sequence in river dikes can be very discontinuous as a result of dike breaches.
It can be inferred from the sparse excavations of dikes ( Figure 1) and from the more numerous cores and CPTs, that many dikes have been constructed during several phases of reinforcement and that often soil material was used that was chosen on the basis of availability rather than its engineering properties. However, their history is often undocumented and observations on dike composition are insufficient to reconstruct a single dike transect, let alone a spatial continuous mapping. Consequently, the stability of the river dike slopes is often very uncertain due to poor estimations of soil properties and groundwater pore pressures. Currently, only average properties of dike material are often derived, or a global statistical uncertainty is applied in a probabilistic analysis (Cao et al., 2016;Gui et al., 2000) to assess its probable response on dike hydrology and stability. Although this approach might be applicable for slope stability problems (van der Krogt, 2022), it ignores the internal structure and heterogeneity of the dike, in which a weak zone can significantly decrease dike macro-stability (Tabarroki et al., 2021). In addition, a wedge of sandy material can increase seepage and pore pressure heads, decreasing the shear strength of the material. Thus, to improve both pore pressure estimations and slope stability calculations internal heterogeneity needs to be included in dike buildup reconstructions.
Two complementary strategies can be deployed to increase information density and improve dike buildup reconstructions. The first strategy is to gather more fieldbased research, using, for example, geological cores and cone penetration tests (CPTs), but the resulting detailed point observations on layer depths and material properties are insufficient for a full reconstruction and in many cases, no local ground truth data of the dike is available at all. The second strategy focuses on point-data interpolation to simulate lithological sequences. Most common are the pixel or voxel based techniques, for example, kriging (Isaaks & Srivastava, 1989), Markov Chain (Elfeki & Dekking, 2005;Li et al., 2004) and multiple-point geostatistics (Feyen & Caers, 2006). Other techniques are object-based methods, which place preconstructed units with certain characteristics in the model domain consecutively (Deutsch & Wang, 1996) or process-based techniques, which simulate the construction or sedimentation process that formed the geological sequence (Pyrcz et al., 2009). However, most of the pixel or voxel-based techniques are less suitable for river dikes simulation given their very diverse and distinct, manmade lithological sequences. Moreover, some techniques depend on abundant training data and local point data for simulation conditioning (Tahmasebi, 2018), that is virtually absent for dikes.
Thus, any proposed method for simulating dike heterogeneity should be capable of simulating the internal structure and composition of dikes and overcome the lack of ground truth data by using general knowledge on dike buildup. In addition, the method should be able to incorporate ground truth data into the simulation, if any is present. These prerequisites require a combination of object-based and process-based techniques, as the process of iterative dike reinforcements can be understood from the general history of river dikes and object-based methods can most easily incorporate information on structural geometry. F I G U R E 1 Example of historical dike buildup as seen in the archeological excavation in Vlaardingen (schematized in Figure 4, see Appendix A for source). The variation in color is a first indication for variable material textures, which are mostly variations of clay-silt mixtures in this dike, with sandy to gravel deposits on the top, which form the road foundation In this study, we present a model to simulate multiple, equally probable versions of a dike's internal stratigraphy. The model uses general information to simulate layer geometry and lithology and can be conditioned on ground truth data where available. The goal of this study is to analyze the performance of this novel method for simulating nonnatural discrete lithological sequences in dikes, which would improve pore pressure estimations and dike stability calculations when compared with methods using average properties with statistical uncertainties. We also apply the method to existing dikes to evaluate the effect of dike interior heterogeneity on groundwater flow under elevated river levels, as a key condition for dike failure. This application highlights the added value of a probabilistic representation of dike buildup in the current practice of dike safety assessments.

| General methods
We developed DETRIS (Dike Erection Tessellation using Regionally Inherited Statistics) to simulate dike buildup, which uses general information on layer geometry and lithology from archeological excavations of dike crosssections (Appendix A; Table A1 and Figure A1). The performance of DETRIS is analyzed by its ability to simulate realistic dike buildups in terms of layer geometry and lithology types. The applicability of DETRIS for dike safety assessment is analyzed by coupling the algorithm to a groundwater flow model and dike safety and comparing these results to an expert assessment of dike slope stability.

| Process based dike construction based on historical assumption
Methods of dike construction in the Netherlands differed substantially throughout history, leading to different heterogeneity patterns of dike buildup. In addition, the inner (landward) and outer (riverward) sides of the dike are known to be constructed of different material: the outside needs to be resistant to wave erosion and therefore often mostly consist of clayey material, whereas the inner side is often a low sloping berm of more sandy material.
For this study, we distinguish four time periods with distinct construction methods (see Introduction) and sufficient sample points to derive meaningful characteristics.
A first period until the closure of the river dike network (1500 CE), a second period until the emergence of largescale dike reinforcements (1800 CE), a third period covers further institutionalizations and the start large scale reinforcement projects and a fourth and final phase comprises modern reinforcement techniques starting from 1950. This last phase is often not present in the sampled dikes as it was of smaller interest from an archeological point of view. We distinguish three geometrical zones along the dike: its outer (river-facing) side, its crest, and the inner side. The dike crest encompasses the top of the dike at those sample points where slopes are below 1:6 (vertical over horizontal distance). Sample locations with an absolute slope value >1:6 belong either to the outer or inner side. Creating geometrical zones based solely on slope is found sufficient, although some errors arose in the automatic classification. However, it has the advantage of rapid pre-processing of new data in case additional crosssections become available. By intersecting the temporal and geometrical zones, a total of nine (3 periods Â 3 dike zones) combinations is created (Figure 2c), each with distinct characteristics.

| Deriving general information of dike characteristics
General statistics of dike buildup characteristics are derived from historical cross-sections of dikes, which are digitized from multiple archeological reports (Appendix A). These cross-sections contain both high resolution layer geometry, detailed material description and dated ages of some layers. At a 10 cm spacing, sample points are created, at which the intersecting layer is selected ( Figure 2). First, the layer material type at that location is determined. Second, the layer thickness (H) is sampled by calculating the length of the shortest line through the sample point connecting both sides of the layer. Third, the layer slope (α) is calculated as perpendicular angle of that line with the horizontal. All sampled points are grouped in the corresponding combination of period and geometry ( Figure 2c). For each combination, the statistics of layer material type is summarized into histograms and the statistics of layer thickness and layer slope are re-sampled to a continuous probability density functions using a kernel density estimate (Scott, 2015).
The statistics of layer material type, layer thickness and layer slope are derived per cross-section, and thus provide local information on dike buildup. Furthermore, regional statistics are derived by aggregating the statistics of all separate dikes. These regional statistics can be used to inform the model even at locations without local statistics on dike buildup.

| DETRIS algorithm
We incorporated the general statistics of dike buildup characteristics per dike construction phase in DETRIS to hierarchically reconstruct the buildup of dikes in an object-based fashion. DETRIS starts with selecting either all local or regional statistics for slope, thickness, and material type. The local statistics are selected if they are available for that dike or a dike in the direct vicinity. Then, a simulation space is created (Figure 3) from the oldest known construction time-period. If no information of construction periods is present, the current dike surface is used. Subsequently, the statistics corresponding to that time-period are selected. An initial dike is created and placed centered at the bottom of the simulation space.
DETRIS iteratively adds a new reinforcement layer on top of the already present dike until the simulation space is entirely filled. In this iteration, first the next reinforcement location (crest, inner, or outer) is selected on the principle of spatial accommodation, that is, the larger the available space between the pre-simulated dike and the final know dike surface, the larger the chance that the next reinforcement stage will be on that location. Then, a new layer is created by sampling a thickness, slope, and material type from the statistics for the reinforcement location and added, until the outline of a dike construction phase is reached. If this is an intermediate phase, the statistics for slope, height and material type are updated to meet those of the next construction phase. If the final outline of the dike is completely filled, the process is completed. In this way DETRIS simulates one possible realization of dike buildup. By creating multiple samples from the same probability density functions, DETRIS can generate many realizations of dike buildup, which are all equally probable.

| Conditioning to ground truth data
If local point observations (e.g., cores, interpreted cone penetration tests) of the dike are available, each newly added layer that intersects with this ground truth data is assigned the observed thickness and material type ( Figure 3). In case the thickness of the simulated layer is smaller than that of the ground truth data layer, the original thickness is kept, otherwise the thickness is reduced to that of the lowest intersecting ground truth data layer. If a simulated layer intersects with multiple ground truth datapoints, it must comply with all.

| Cross-section descriptions
Two cases are selected to analyze the performance of the DETRIS algorithm, and one case is selected to explore the added value in dike slope stability assessments.

| Cross-sections for analysis of DETRIS performance
To analyze the material and shape heterogeneity and uncertainty related to DETRIS, two known cross-sections from the database (Appendix A) are selected: Vlaardingen (Maassluissedijk Vettenoord) and Lent (Bemmelsedijk). These are amongst the most complete historical crosssections available and are located, respectively, in the more distal and proximal parts of the Rhine delta ( Figure 4). The Vlaardingen dike is therefore mostly constructed of clay material, while the Lent dike has a strong variation of clays, sands, and silts. Furthermore, the Lent dike was almost entirely constructed before 1500 CE, while the Vlaardingen dike was raised and widened considerably after that time.
F I G U R E 2 Derivation of dike statistics from archeological dike cross-sections. First, a grid of sample points is created (a) at which the layer material type, layer thickness (H) and layer slope (α) is selected (b). These statistics are grouped given their age and location in the dike. For example, one group consists of points that were constructed as dike crest (upward triangle) between 1500 and 1800 CE (blue), as denoted by the combined symbol (blue upward triangle) in C. Here, the outer slope refers to the side of the dike next to the river, the inner slope is that on the side of the protected area (polder).

| Cross-section for DETRIS application
The application of DETRIS in relation to the assessment of groundwater conditions and dike slope stability is shown for the northern Lek dike east of the village Nieuwegein (Figure 5a,b). A survey resulted in a complete lithological characterization of the natural subsurface ( Figure 5c) and some point observations (ground truth data) on the dike interior from a cone penetration test (CPT, Figure 5d).

| Performance of DETRIS
The performance of DETRIS is analyzed for the Vlaardingen and Lent dikes (Figure 4). The performance is assessed in terms of the ability of DETRIS to create similar patterns of heterogeneity as found in historical dikes, and the ability of DETRIS to simulate the correct material texture given variable quantities of available ground truth data.
In the assessments, local statistics of slope, layer thickness and material types, that is, from the corresponding cross-section, and regional statistics, that is, the aggregated over all cross-sections in the database, were evaluated. For each dike and selected set of statistics DETRIS created 100 hypothetical dike realizations using the dated outlines of the cross-sections for the construction phases.

| Heterogeneity measures
The performance of DETRIS related to material and layer shape heterogeneity is analyzed on four simple quantitative heterogeneity metrics: fractal dimension D (O'Neill et al., 1988), relative contagion RC (Li & Reynolds, 1993), relative evenness RE and relative patchiness RP (Romme, 1982). They are frequently used in landscape characterization (Li & Reynolds, 1993) and represent different components of spatial heterogeneity. Fractal dimension (D) measures shape irregularity and is calculated using: where A k and P k are the area of perimeter of layer k. D is derived from a linear regression between the natural logarithms of A k and P k over all patches. Higher values indicate an increase in patch shape irregularity.
The relative contagion (RC) measures patch continuity and spatial arrangement and is expressed as: where n is the number of different lithologies present and P ij is the probability that two adjacent pixels consist of lithology i and j. Higher values indicate a larger clustering of same type patches or a larger patch size. The relative evenness (RE) specifies the proportions of each lithology type and is expressed as: where P i is the probability of a random pixel to be of type i. Higher values indicate a more uniform patch type distribution.
The relative patchiness (RP) is a measure of the difference between adjacent layers and is expressed as: where E ij is the number of transitions from type i to type j, T is a unit chosen to indicate the dissimilarity for patch types i and j and N is the total number of transitions. The F I G U R E 4 Overview of locations and cross-sections of the two known historical dikes used to measure the performance of the DETRIS algorithm. The age lines correspond to the location of the dike surface at that time (years CE). The presented heights are relative to the Dutch ordnance datum. Material types according to USDA classification system (Soil Science Division Staff, 2017) dissimilarity value T equals log 10 K sat ð Þ corresponding to the material lithology (Table 1), as it is assumed representative for material type differences. Higher values indicate more and larger transitions between patches.
These four metrics are compared between the simulated and observed cross-sections, and a combined performance is measured by the normalized root mean squared error (NRMSE) by: with N as the number of heterogeneity measures (four, as described above), y n the average simulated value of a specific heterogeneity measure and RMSE n is the root mean squared error of that heterogeneity measure, calculated as: with T the total number of model runs, b y n the observed value for the heterogeneity measure and y t,n the simulated value for the heterogeneity measure.

| Uncertainty quantification
The balanced accuracy (BA) is used to quantify the uncertainty of the DETRIS model. It can be used to address spatial differences (i.e., within the simulated dike) in accuracy, but also to address differences between multiple runs with different forcing. If each class is equally well predicted, or if there only is a single true value, this term is equal to normal accuracy. Otherwise, the value decreases as each sample is weighted according to its inverse prevalence. It is written as in which TP equals the count of true positives, FN equals the count of false negatives, etc.

| Application of DETRIS for dike stability assessment
The application of DETRIS is evaluated on a crosssection for which a dike slope stability assessment is performed, including estimates on groundwater conditions. As described in Section 2.3.2, this case study is located at Nieuwegein ( Figure 5) and contains ground truth data (a CPT) to which DETRIS will be conditioned.

| CPT classification
The CPT is classified to five lithological classes based on Been and Jefferies (1993): sands, sand mixtures, silt F I G U R E 5 Overview of case study location and cross-section. The maps on the top row show the case location within the Netherlands (a) and a more detailed overview of the cross-section (b). The bottom row (c) shows the cross-section at that location, given the surface elevation and subsurface properties (Table 1). A zoom of the dike (d) is presented that focuses on the outline of the dike and the classified cone penetration test (CPT), which is used as ground truth data by the DETRIS-algorithm. The presented heights are relative to the Dutch ordnance datum van WOERKOM ET AL. mixtures, clays and organic soils (Table 1). These material types are used to simulate the dike buildup using DETRIS, and the originally used material type histogram (Figure 7) is reclassified to these units.

| Basis of expert scenarios of groundwater conditions and dike slope stability
The available expert dike slope stability assessment is performed with the Dutch guidelines (Rijkswaterstaat, 2021) for dike slope stability assessment. These guidelines of this standard (TAW, 2004) specify two scenarios for the phreatic groundwater table in the dike and one for the pressure head in the underlying sandy aquifer.
The expert groundwater scenarios are imported in slope stability software D-Stability (Deltares, 2019). The Uplift-Van limit equilibrium method is used for slope stability calculations, which assumes a dual circular slip plane. The slip plane consists of an active circle, a passive circle and a compression bar in between (Van, 2001). The slip plane has to start in the crest or the upper half of the landside slope. A C-Phi model for calculating shear stresses is used, and each layer in the DETRIS simulations is assigned a characteristic value of cohesion (c'), friction angle (φ), material weight (γ) and saturated material weight (γ sat ) given its simulated material type. The expert scenario dike is assumed homogeneous and is assigned a single value for each of the parameters (Table 1).

| Calculating groundwater conditions and dike slope stability from DETRIS
The selected Nieuwegein cross-section ( Figure 5) has no archeological observations (Appendix A) in its vicinity, so the DETRIS realizations are based on the regional statistics. The classified CPT serves as ground truth data. DET-RIS generates an ensemble of 100 possible dike realizations. For each realization steady-state pressure heads are calculated, and a dike slope stability calculation is performed on the DETRIS simulation and the resulting pore pressure heads.
The steady-state hydrological simulations of the entire cross-section are performed using the MODFLOW 6 software (Langevin et al., 2019). The 2D hydrological model is setup using a cell size of 0.2 m both vertically as horizontally and includes the dike and the subsurface (Figure 5c). At the left end the river has a connection with the lower aquifer. Here the imposed river stage constitutes a head-controlled boundary condition (Figure 2) handled by the MODFLOW river package (Hughes et al., 2017). On the inner side of the dike head-controlled conditions ( Figure 6) are handled by the MODFLOW drain package that creates seepage points that permit T A B L E 1 Subsurface types used for case-study and corresponding hydrological and geotechnical parameters Note: Values of saturated conductivity are derived from (Pachepsky & Park, 2015) and geotechnical parameters are based on EC7 (CEN, 2004).

Material type
F I G U R E 6 MODFLOW cross-section and boundary condition types. Boundaries without specified boundary conditions are no-flow boundaries. The different shades of gray are based on the lithology, which is more explicitly visualized in Figure 5c outflow only (Hughes et al., 2017). A drainage ditch with a constant head 0.5 m below the surface is located on the landward boundary. The subsurface material types and simulated dike material types are assigned a saturated conductivity (K sat ) based on characteristic values (Table 1). In the steady-state simulation, the river is at a high stage, with a head of 6.73 m with respect to ordnance datum (recurrence interval of 1:30,000 years, the safety standard for this dike stretch). The dike slope stability calculations are again performed in D-Stability software (Deltares, 2019). The DET-RIS realization replaces the expert schematization of the dike material, while the natural subsurface is kept the same ( Figure 5). Each layer in the DETRIS realization is assigned the C-Phi parameter values corresponding to its material type (Table 1). In addition, the resulting pressure heads from the steady-state MODFLOW model are inserted in the D-Stability software. From the DETRIS realization and corresponding pressure heads, the inner (landward) slope stability is expressed by the safety factor F (see supplementary material Data S1 calculation for example).

| General characteristics of archeological dike cross-sections
Out of the 66 digitized cross-sections, 12 cross-sections had sufficient detail to derive general statistics of dike buildup characteristics in the form of probability density functions for layer slope, layer thickness and layer material type. The layer slope is centered around a slope of zero, with negative values representing inner side slopes and positive values representing outer side slopes (Figure 7). The median layer slope is exactly zero, indicating inner and outer slopes are generally equally steep.
The absolute median layer slope equals 0.21, roughly being equal to a vertical to horizontal ratio of 1 to 5. The median layer thickness is 0.50 meter, with 5 and 95 percentiles at 0.13 and 1.94 m. In our digitized cross-sections most dike layers consist of clay and loam (71%), with the largest remainder consisting of sand or gravel (26%). However, peat, wood and stone are also found in the historical dikes (3%).

| Visual DETRIS example
To illustrate how DETRIS simulates construction history by consecutively selecting new reinforcement layers, the dated construction of the real Lent dike (Figure 8) is presented next to an arbitrary DETRIS realization of the same dike. The total construction period is divided in two phases: a first phase with construction and reinforcements made up to 1500 CE and a second phase with reinforcements made after 1500 CE. As outlined, these construction phases are linked to the different statistics to parameterize DETRIS, resulting in two separate simulation spaces that can be traced across the realizations.
The early reinforcements particularly consist of clayey to silty material in both the real dike and the DETRIS dike. The real dike equally extends the reinforcements to either side of the dike center, confirming our DETRIS assumption that the position of any next reinforcement is based on the principle of spatial accommodation, that is, that there is no clear preference to one side. For the second (youngest) selected temporal period, both the DET-RIS dike and the real dike have a larger occurrence of sand and gravel reinforcement layers. The DETRIS realization visually seems to mimic the construction history and resulting heterogeneity of the real dike relatively well, but a statistical analysis is needed to demonstrate this similarity on a larger set of realizations.
F I G U R E 7 Histograms of regional dike properties from all high-detail cross-sections (Appendix A). The sampled properties are layer slope, layer thickness and layer material type. Negative layer slopes indicate inner side slopes and positive layer slopes indicate outer side slopes. Material types according to USDA classification system (Soil Science Division Staff, 2017)

| Performance of DETRIS algorithm
The performance of DETRIS in probabilistic terms is assessed on two known cross-sections (Lent and Vlaardingen, Figure 4), using statistics from both local as regional statistics. The performance is based on two indicators: The ability of DETRIS to create similar patterns of heterogeneity as found in historical dikes, and the ability of DETRIS to simulate the correct material texture given variable quantities of available ground truth data.

| Performance of simulating dike heterogeneity
Four simple quantitative heterogeneity metrics (Section 2.4.1), being fractal dimension (D), relative contagion (RC), relative evenness (RE) and relative patchiness (RP), are compared between 100 DETRIS realizations and the two known dike cross-section. The average normalized root mean squared error (NRMSE) of all heterogeneity indices between all the realizations and real the values are seen as a general indication of model performance (Figure 9).
The NRMSE using local and regional statistics for dike buildup characteristics, respectively, are 0.18 and 0.21 for the Lent case and 0.20 and 0.23 for the Vlaardingen case. This indicates that a dike can be more precisely reconstructed by DETRIS if the local statistics of that dike are known. For both Lent and Vlaardingen cases, DETRIS performs best in the fractal dimension (D) measure, indicating that the simulated layer shapes match well with the layer shapes in the real dikes. For the tested cases, the RP of the local statistics generally tends to be low and of the regional statistics high, indicating that layer transitions given the regional statistics tend to be associated with shifts in material properties that are more abrupt than observed in reality. In line with this the Vlaardingen RE is often too high, which might suggest that DETRIS is less suitable for simulating more homogeneous cross-sections.

| Performance of model accuracy and effect of ground truth data conditioning
The DETRIS-algorithm can also be conditioned on ground truth data, often available as vertical cores. These cores are often placed at the outer toe, at the dike crest or at the inner toe of the dike. Three scenarios are compared here: a scenario without ground truth data (Figure 10a,d), one location of ground truth data ( Figure 10e) and three locations of ground truth data (Figure 10c,f) for DETRIS conditioning. For each scenario, 100 DETRIS realizations are performed using local and regional statistics for both Lent and Vlaardingen. The spatial balanced accuracy (BA) is only assessed on 100 realizations of Lent with local statistics (Figure 10d-f).
Independent of the use of local or regional statistics, the median BA increases when adding ground truth data, while the interquartile range of BA decreases. The local statistics result in a higher BA than the regional statistics, but this difference decreases when adding one core and almost diminishes when three ground truth data cores are used. There is a large difference in BA between the two presented cases, Lent and Vlaardingen, where the more uniform Vlaardingen case shows the largest relative improvement from no cores to one core, while the more heterogeneous Lent case shows the largest relative improvement from one core to three cores.
The BA increase of the DETRIS model is mostly caused by an increase in the proximity of the available ground truth data (Figure 10d-f), with BA > 0.75 for on average 1.7 m on each side of the ground truth data. At locations further away from the ground truth data the BA stays the same or can even decrease compared with the noninformed (no cores) reconstruction.
F I G U R E 8 Historical evolution of a real dike (Lent) based on the archeological excavation, and a realization of dike construction history by the DETRIS algorithm using the local statistics from the observed dike. Total construction period is approximately from 1200 to 1950 CE. The construction history is divided in two phases, separated by the dotted line, representing the historical dike surface at approximately 1500 and 1950 CE

| Implications of DETRIS for groundwater levels and dike stability
The application of DETRIS is shown for a different crosssection (Nieuwegein, Figure 5). At this location, the regional statistics are used as no local information is available except a single CPT. For this location, two expert scenarios on groundwater conditions and the corresponding safety factors (F) for dike slope stability are available (see Section 2.5.2).
Based on the 100 dike buildup realizations by DET-RIS, the simulated pressure heads in the aquifer beneath the dike are in general agreement with the schematized expert scenario ( Figure 11B). For the realizations, the variation of dike buildup material type results in a high variation of the simulated water tables in the dike. At the outer crest location (width = À3.5) the difference between the 5 and 95 percentile of the simulated phreatic lines equals 1.59 m, while the two expert schematized phreatic lines differ 0.98 m. Around the CPT location and at the inner toe of the dike, the variability decreases. At the inner crest location (width = 3.5) the difference equals 0.48 m for the simulations and 1.05 m for the schematized phreatic lines, and this decreases even further towards the inner toe, as the dike is drained by the sandy wedge present. At the CPT location, the reduced F I G U R E 9 Comparison of heterogeneity indices of DETRIS algorithm based on local and regional statistics of dike buildup characteristics and the heterogeneity indices of the actual dikes: Fractal dimension (D), relative contagion (RC), relative evenness (RE), and relative patchiness (RP) F I G U R E 1 0 Performance of model predictions given one or multiple ground truth data locations for conditioning. The top row (a-c) shows the balanced accuracy (BA) on both cross-sections and local/regional statistics, the bottom row shows the balanced accuracy (BA) spatially of the Lent dike based on local statistics 5-95 percentile range of groundwater pressure heads emphasizes that the effect of any ground truth data extends beyond decreasing the uncertainty in material texture (as shown in the previous section, Figure 10).
The two safety factors F for dike slope stability on the inner (landward, right side) slope are 1.12 with a high phreatic level and 1.53 with a lower phreatic level (Figure 11c). Based on the DETRIS realizations, the median dike stability F is 1.41, with a 5 and 95 percentile of 1.35 and 1.51. Thus, the simulated safety factors F are in between the two expert schematizations and have a smaller variation. Overall, the DETRIS-based values are nearer to the safety factor resulting from the schematized groundwater scenario with a lower phreatic level as the draining conditions near the inner crest strongly control the phreatic levels in this case.

| Relevance of local statistics for dike heterogeneity assessments
To assess the dike safety on a longer dike stretch, information on dike material variation along the dike on a small scale (several meters) is crucial. However, to our knowledge only few observations exist from which the parameters needed for DETRIS can be derived, and these F I G U R E 1 1 Results and implications of using DETRIS for dike stability simulation. The entire cross-section is shown in Figure 5. For 100 dike realizations using DETRIS (of which four realizations are shown as example [panel a]), the groundwater level (b) and dike stability (c) is calculated and compared against schematized expert scenarios. The simulated phreatic levels (b) (light blue) are mostly in-between the two expert schematized phreatic levels (dark blue). Similarly, the resulting dike stability safety factors (F, panel c) are also mostly in between the safety factors resulting from the two expert schematized phreatic levels (dotted lines) are often far apart. However, an indication of the spatial correlation of material properties is essential as it indicates how representative the local and regional statistics are that are used to parameterize DETRIS.
The spatial variability of dike buildup is analyzed based on the material types of all 66 cross-sections in the larger database (see Section 3.1), which include both the 12 detailed profiles used for the parametrization of DET-RIS as well as other less detailed reports on historical dikes. Two parameters are assessed: the previously discussed RE, as a measure of dike material heterogeneity, and the average saturated conductivity K avg based on Table 1, as a measure of average dike properties. The spatial similarity of these metrics between dikes is assessed using variograms (Matheron, 1963). A spherical variogram model is fitted including the nugget effect, with a maximum lag distance of 100 km and a uniform binning strategy.
The K avg has a nugget of 0.989, which is 50.5% of the total sill of 1.958, and an effective range of 63.4 km. The RE has a nugget 0.033, which is 53.3% of the total sill of 0.063, and an effective range of 7.7 km ( Figure 12).
The variograms of Figure 12 suggest that the use of local dike statistics is only useful within 7.7 km from a known dike cross-section, as at further distances no spatial similarity in dike buildup may be expected. Of the approximately 1850 km of primary river dike in the Netherlands, 17% is located within a radius of 7.7 km from a detailed dike cross-section. As the local statistics from detailed dike cross-sections substantially decrease the uncertainty in dike buildup and groundwater conditions, increasing the number of detailed cross-sections from which local statistics can be derived, would increase the overall accuracy. However, the high nugget-to-sill ratios ($50%) indicate that even on smaller scales large variations can occur, and that ground truth data, in combination with a successfully applied interpolation or simulation model, is also needed for accurately predicting dike buildup. This was also shown by the balanced accuracy (BA), where for Vlaardingen, local statistics without ground truth data lead to an median balanced accuracy (BA) of 0.17, while regional statistics with three ground truth data cores for conditioning lead to a BA of 0.49 and local statistics with three ground truth data cores lead to a BA of 0.51 ( Figure 10). Thus, the model benefits more from adding ground truth data at multiple locations when regional statistics are used, than it benefits from adding a limited number of cross-sections to provide the local parameterization. Also, cross-sections are very destructive and often not allowed in primary dikes. However, we recommend a more intense collaboration with archeologists to increase the density of high-quality data on dike material buildup to inform dike reconstructions, in particular as the delineation of the reconstruction phases is an open question in our method. Again, replacing it by aggregated statistics for the entire dike buildup and using ground truth data can solve this partly, but the results will deteriorate considerably without ground truth data conditioning.

| Dike safety assessment improvements by DETRIS
DETRIS can provide more realistic heterogeneous dikes than those currently in use, as it does not assume uniform material characteristics in the entire dike. Its fast and effective algorithm can provide multiple realizations of dike buildup, based solely on input statistics of material and layer characteristics, and optionally incorporate ground truth data. As such, it produces a data-based approximation of the uncertainty in dike failure risk as a function of dike heterogeneity and groundwater conditions. This probabilistic approach results in a better informed assessment of dike slope stability. In addition, it F I G U R E 1 2 Variogram models of two dike properties. Both show a large nugget, indicating local variability at scales smaller than sampled in this research is significant can easily incorporate (more) local statistics, ground truth data, or knowledge on dike construction history in an objective manner, contrary to human-controlled expert schematizations. In this case, the probabilistic assessment resulted in a smaller uncertainty when compared with expert judgment (Figure 11). While this may vary from case to case, it reflects the strength of DETRIS in evaluating the effect of the buildup and the uncertainty therein on the stability of dikes in a mechanistic manner.

| Suggestions for model improvements
DETRIS is currently capable of simulating realistic dike buildup on a 2D cross-sections, but the algorithm is prone to simulate too large differences in composition between adjacent layers. In the observed cross-sections similar layers were often more clustered. Like conditioning on g point data, adding conditional probabilities to represent inter-layer dependencies can result in more realistic dike buildup realizations. The use of recent developments in continuous geophysical data interpretation, such as electrical resistivity tomography (Chavez Olalla et al., 2021) or ground penetrating radar (Di Prinzio et al., 2010) are promising to increase the availability of ground truth data too. In addition, coupling our object-based model to other geostatistical methods, such as multiple-point geostatistics (Tahmasebi, 2018) can provide more detailed indications of material uncertainty across but also along the dike. Such a coupling has already been successfully applied to simulate natural deposits (Michael et al., 2010).

| Outlook and suggestions for further research
The explicit characterization of the legacy of dike construction on the material composition as implemented in the DETRIS algorithm can assist in obtaining better informed safety calculations and improve our understanding of the governing processes on different temporal scales. An example of a short-term process as presented in this paper is the location of the phreatic line during high river water levels and the dike slope stability. Other stability related features, such as the presence of an erosion resistant core or small sand lenses in the dike can be better assessed using DETRIS as well. They are important as a source of residual strength after dike slope failure and for having a major influence on the intrusion length of pore pressures.
In terms of long-term processes such as soil formation and weathering, DETRIS can link the composition of a particular layer and its properties to material and age. In addition, compaction of the layers can change the dike buildup geometry and material properties. This, eventually, may alter the behavior of the dike during prolonged periods of droughts due to irreversible shrinkage of clay and peat and the formation of fissures, which could affect its ability to withstand high water levels that follow, thus introducing hysteresis.

| CONCLUSION
In this work, we developed the Dike Erection Tessellation using Regionally Inherited Statistics (DETRIS) model to simulate dike buildup based on historical dike cross-sections. By applying this model, a new valuable source of information is added to improve estimates of dike material and its heterogeneity, and provide uncertainty estimates of flood safety of the corresponding dikes.
• The implementation of a process-based model (DETRIS) enables a probabilistic assessment of dike interior buildup based on observed characteristics of historic dikes, which goes beyond the current assumption of homogenous dikes. • Dikes have a high internal heterogeneity which is captured by the model. Local information improves the predicting capacity of dike buildup, but regional statistics, nonetheless lead to representative dikes. • Dike cross-sections only have a predictive capacity within several kilometers, and considerable variation is observed even on smaller scales. Conditioning to ground truth data cores or CTPs provides a good alternative at locations outside this range. • Using DETRIS-based scenarios of dike buildup to simulate groundwater levels and dike stability provides a more comprehensive way of incorporating heterogeneity than arbitrary expert-based estimations and narrows down the uncertainty of dike slope stability analyses while incorporating local and adverse effects of more permeable or less competent layers. • DETRIS is capable of simulating nonnatural discrete lithology sequences, and it can be further improved by adding layer inter-dependency or possibly coupling it to a multiple-point geostatistics algorithm. • We recommend a more intense collaboration with archeologists and responsible institutions for dike maintenance to increase the density of high-quality data on dike buildup material and ages of dike erection phases than currently present.

ACKNOWLEDGMENTS
We greatly appreciated the support of Martin van der Meer and Johan Hockx providing data on the case-study and ideas on the application of the algorithm. We also greatly thank the late Michel Lascaris for providing references on archaeological dike cross-sections.

DATA AVAILABILITY STATEMENT
Data sharing is not applicable to this article as no new data were created or analyzed in this study.
ORCID Teun van Woerkom https://orcid.org/0000-0001-7629-3549 F I G U R E A 1 Location and level of detail of archeological sources of dikes in the Netherlands. Those with sufficient quality are listed in the table  above