Geophysical tomography as a tool to estimate the geometry of soil layers: relevance for the reliability assessment of dikes

ABSTRACT The geometric variability of soil layers is a large source of uncertainty in the reliability assessment of dikes. Because direct samples of the subsurface soils are often insufficient to capture the complexity of the subsurface, geophysical methods provide a powerful source of complementary information. A combined approach to estimate the geometry of soil layers is presented. The approach combines local point data, i.e. data obtained from a CPT or a borehole log, and geophysical tomography in a universal cokriging framework. The approach uses the contact points between soil layers obtained from local point data and the orientations of the layers derived from geophysical tomography. To reduce subjectivity in the interpretation of tomographic images, an automated edge detection technique was used. The combined approach was applied to characterise two test sites where the presence of paleochannels locally change the geometry of soil layers. The results show that a combined approach enables the reduction of sampling efforts with an improved estimation of geometric variability.


Introduction
Dikes form an essential part of the primary flood defences along the coast and major rivers in the Netherlands. Multiple failure mechanisms threaten the stability of these dikes. Hence, they are subjected to periodic reliability assessments (MinIM, 2016;De Waal, 2018). In reliability assessments, a special focus lies on the geological schematisation of the subsurface (Hijma and Lam, 2015). Characterizing geometric variability of soil layers is a key step within the geological schematisation process. In the failure mode of macrostability, for example, the thickness of a weak layer determines the shape of the slip surface and the reliability of the dike. Geometric variability is also important for the reliability assessment of dikes in terms of piping. In clay-over-sand dikes, the thickness of the clay layer in the hinterland generates resistance against uplift which is the first phase of piping (Sellmeijer and Koenders, 1991).
The Dutch subsurface is notoriously heterogeneous. It is built up by sequences of gravely sandy aquifers alternated by confining clayey aquitards. The uppermost aquifer, which was deposited during the last glacial period (Pleistocene epoch), forms a nearly continuous sandy substrate. These deposits are overlain by the Holocene deltaic wedge which forms a heterogeneous confining layer. This layer is mostly composed of aquitard floodplain clays, clay-fine sand, dominated intertidal flats, and peats. However, it is dissected by (partly) isolated alluvial and tidal channel sand bodies (Bierkens, 1994;Weerts, 1996;Hijma and Cohen, 2011). These channel-belt sand bodies act as shallow-depth aquifers and locally occur directly underneath dikes. In addition, small-scale variability of sand body architecture also occurs which is the result of autogenic processes. Due to this complexity, it is challenging to characterise the geometric variability of soil layers with conventional site investigation methods. Methods, such as the Cone Penetration Test (CPT) and borehole drilling, sample the subsurface in detail, yet locally. In the Netherlands, CPT samples are normally collected with a spacing of 100 m (ENW, 2012). In the Dutch subsurface; however, the geometry of soil layers often varies on scales smaller than 25 m (Hijma and Lam, 2015). The mismatch between site investigation density and actual subsurface variability leads to schematisation uncertainties. Although schematisation uncertainties are compensated with safety factors (ENW, 2012), large uncertainties result in uneconomical designs or unnecessary reinforcement of existing structures. Therefore, it is of paramount importance to reduce these uncertainties without exhaustively sampling the subsurface. For this purpose, geophysical methods are a well-established option because they map the subsurface in a horizontally-continuous manner.
Even though geophysical methods provide valuable insights into the variability of soil layers, their use has been limited to the visual interpretation of geophysical images. Few attempts have been made to objectively derive information from geophysical data. For example, deterministic (Auken and Christiansen, 2004) and probabilistic (de Pasquale et al., 2019) inversion methods have been used to estimate the geometry of soil layers and material properties. Alternatively, Hsu et al. (2010) and Chambers et al. (2012) applied automated detection techniques to tomographic images in order to estimate the geometry of soil layers. The problem is that geophysical data are often affected by instrumental drift, lateral heterogeneity, and lack of resolution (Minsley et al., 2012;Delefortrie et al., 2019). In such cases, the quantities estimated with these methods, e.g. geometry, though informative of the variability trend, are inherently inaccurate.
We propose a combined approach to estimate the geometry of soil layers. We combine the trend of geometric variability estimated from geophysical data with accurate data derived from boreholes. We use the potential field method (Lajaunie, Courrioux, and Manuel, 1997) to combine both data sets. The potential field method generates a geometric model of the subsurface via universal cokriging of the layer orientations and contact points between soil layers. We use the Laplacian edge detection technique to estimate soil interfaces from geophysical tomography. The layer orientations are then calculated as the dip angle of the tomographic interfaces. Thus, we use layer orientations derived from tomography and the contact points between soil layers obtained from directly sampling the subsurface. We test the approach with two electromagnetic geophysical methods that are widely used in the site investigation of dikes, namely Electric Resistance Tomography (ERT) and Electromagnetic Induction (EMI).
We first describe the proposed approach and the main steps behind it. We describe the construction of geophysical tomography, the technique to obtain layer orientations from tomography, and the potential field method. Afterwards, we present a proof of concept of our approach with a synthetic two-layer model with internal variability. Our approach is then applied to two study sites. One site located in an alluvial environment across an old river channel and the other site located in a tidal environment along the longitudinal section of a dike. In both cases, the geometry of the upper layer is characterised.

Methodology
We propose an approach to estimate the geometry of the soil layers of a geological setting (Figure 1). Local point data and geophysical tomography are used as input in this approach. Local point data, such as CPTs and borehole logs, are used to derive the contact point or interface between soil layers. ERT and EMI tomography (Section 2.1) were used to derive the orientations of the soil layers. The orientations of the soil layers were calculated from the edges automatically detected in the tomographic images (Section 2.2). Finally, orientations and contact points were used to estimate the geometry of soil layers via the potential field method of Lajaunie, Courrioux, and Manuel (1997) (Section 2.3).

Geophysical tomography
Two geophysical methods were considered, ERT and EMI. Although the working principle of each method is different, both methods map the electrical resistivity of the ground. ERT is based on electric conduction while EMI is based on electromagnetic induction. The ERT method allows for fine-tuning the sensor positions and acquisition array. Therefore, the depth of investigation and sensitivity of the sensors are adaptable to different applications. Nevertheless, setting up an ERT system is labour intensive and the data acquisition is time-consuming. On the other hand, data acquisition with an EMI system is up to orders of magnitude faster (Binley et al., 2015). The EMI device used in this study, DUALEM-421, has a fixed geometric configuration. Therefore, fine-tuning the depth of investigation and sensitivity was not possible. A detailed explanation of the working principles of both geophysical methods is presented in Appendices A.1 and A.2, respectively. A tomography, which is a model of the subsurface, was constructed with geophysical measurements. The tomography, in this case, represents the spatial distribution of electrical resistivity in the subsurface. Since tomographic inversion of ERT and EMI is an ill-posed problem, it requires regularisation. Regularisation is an assumption about the form of the tomography. Smoothness regularisation was used to construct the ERT and EMI tomography. A further explanation of the tomographic inversion process is presented in Appendix A.3. The orientations of the soil layers were then estimated from the geophysical tomography (Figure 2(a)).

Layer orientations
The Laplacian edge detection technique was used to automatically detect edges in the tomographic images. The edges in the tomography were interpreted as the interfaces between soil layers. In the Laplacian edge detection technique, edges are defined at the zero-crossings of the Laplace operator applied to the tomography where m(x, y) are the material properties of the medium represented in the tomography i.e. electrical resistivity. Due to tomographic artefacts and geological complexity, fake edges and edges that are not of interest are also detected in the tomography. Thus, additional processing and interpretation are often needed. Processing with smoothing and thresholding is effective in reducing the detection of fake edges (Mlsna and Rodríguez, 2009). Smoothing is effective in reducing the detection of fake edges caused by the small-scale components in the tomography. An appropriate filter strength should be picked based on the frequency content of the tomographic images. On the other hand, thresholding removes edges that show a low tomographic gradient. From the remaining edges in the tomography, the edges that correspond to the geological feature of interest are picked based on the geological interpretation of the tomography. The final selection of tomographic edges is assumed to represent the trend of geometric Figure 1. Schematic of the approach proposed to estimate the geometry of soil layers of a geological setting. The approach uses local point data and geophysical tomography. Local point data are used to derive the contact points between soil layers. Meanwhile, geophysical tomography is used to derive the orientations of the soil layers. The orientations are calculated from the edges detected in the geophysical tomography. The estimated geometry of the soil layers is obtained via universal cokriging of the contact points and the orientations. variability of the soil interfaces. We represent this trend with the orientations of the edges (Figure 2(a)). The orientation is a unitary vector that is perpendicular to the soil layer. The orientations are calculated from the dip angle of the tomographic edges. In two dimensions the orientations become where the dip angle is calculated directly from the tomographic edges where x int and z int are the coordinates of the edges. The dip angle calculated from the tomographic edges is signed which defines the pointing direction of the orientation vector. Alternatively, the orientation vector can be defined in terms of dip and azimuth. The performance of the edge detector depends on tomographic resolution. Resolution in turn depends on measurement physics, regularisation, acquisition design, and physical contrast between soil layers. Appendix 2 describes the performance of the edge detector for several cases of acquisition design and physical contrast in a two-layer model. Appendix 2 shows that lack of resolution leads to poorly mapped regions of the subsurface where automated detection techniques do not perform well. In the regions where automatically detected edges are not informative, assumptions have to be in terms of the contact between soil layers or their orientations. The assumptions made in this part of the process have a large effect on the estimated geometric variability. The advantage of the potential field method (Section 2.3) is that these assumptions are made explicit. Therefore, uncertainties in these assumptions are quantifiable.

Potential field method: universal cokriging
The potential field method (Lajaunie, Courrioux, and Manuel, 1997) was used to combine tomographic orientations with contact points to estimate the geometry of soil layers. The method is based on universal cokriging and offers two main advantages, namely flexibility and objectivity. The method is flexible because it allows two non-collocated variables as input. In the present case, the contact points between soil layers are accurately known from CPTs and boreholes, but they are sparse. Meanwhile, the orientations of the soil layers are known in a horizontally-continuous manner from geophysical tomography, but their location does not coincide with the contact points. The potential field method is also objective because it takes into account the input data sets explicitly. Consequently, the effect of the data and the covariance model on the estimated geometry of the soil layers is quantifiable.
The method generates a potential field which is designated with the random function Z 1 . The potential field represents a proxy to the time of formation. By this definition, the interface between two soil layers was formed at the same time. Also, the gradient of the potential field, Z 2 =∇Z 1 , coincides with the orientations of the soil layers. The gradient points towards younger formations or later formation times. To formulate the cokriging system, the value of Z 1 and Z 2 needs to be known at certain locations. The gradient of the potential field, Z 2 , is known from geophysical tomography (arrows in Figure 2(b)). However, the value of the potential field is not known at the contact points (crosses in Figure 2(b)). Thus, it is convenient to replace the random function Z 1 , which is not known at the contact points, by a new, known, random function defined as where x 0 is a reference contact point which belongs to the same interface as x. To illustrate, in Figure 2(b), x 0 is the first contact point (leftmost red cross) so Z 1new (x) is evaluated at the remaining two contact points. The choice of reference point has no influence in the cokriging estimation of the potential field Z * 1 (x). The function Z 1new equals zero because the value of the potential field along an iso-surface is constant. The potential field method allows for any number of soil layers, but at least two contact points must be known per layer so that Z 1new is defined. The potential field is then generated via universal cokriging of Z 1new and Z 2 . Because of the dependency between Z 1 and Z 2 , the covariance and cross-covariance matrices are derived from a single covariance model, i.e. the covariance model of the potential field. The same is true for the drift functions. Following De La Varga, Schaaf, and Wellmann (2019), the cokriging system becomes where C Z 1new ,Z 1new , C Z 2 ,Z 2 are covariance matrices, C Z 1new ,Z 2 , C Z 2 ,Z 1new are cross-covariance matrices, U Z1new and U Z 2 are drift functions, λ, μ are the cokriging weights, c are vectors which contain covariances and cross-covariances between the existing data points and the interpolation point, and f Z 1new is the vector which contains the drift function of Z 1new evaluated at the interpolation point. The drift function and the relation between covariance and cross-covariance models are elaborated in Appendix A.4. The potential field is estimated at any interpolation point in the domain from the cokriging weights where M is the number of contact points minus the number of soil layers and N is the number of gradients which is a multiple of three. Although the contribution of Z 1new in equation (6) is zero, Z 1new contributes to the cokriging weights associated with Z 2 . The covariance models, which are needed to calculate the covariance matrices and vector in equation (5), can be derived experimentally or heuristically. Experimentally, the covariance models need to be derived from the orientation data (Aug, 2004;Chiles et al., 2004), for the potential field is a mathematical construction not known at the contact points. Since the relation between covariance models is known (Appendix A.4), the covariance models follow from the covariance model of the orientation data. Alternatively, a heuristic approach is often used to define the covariance models. De La Varga, Schaaf, and Wellmann (2019) assume a spherical covariance model for the potential field Z 1 where the variance and range of Z 1 define all the covariance and cross-covariance models. Default values for the variance and range of the covariance function are calculated based on the size of the model domain (Appendix A.4). The estimation variance has no physical meaning when the covariance models are defined heuristically. In that case, the cokriging system can be solved in its dual form which improves computation efficiency (Goovaerts, 1997). The python package GemPy (De La Varga, Schaaf, and Wellmann, 2019) was used to apply the potential field method to the data presented.

Synthetic study site and data simulation
The approach of Figure 1 was applied to a synthetic twolayer model with internal variability. The contact points between soil layers were obtained from sampling the subsurface at two locations. The orientations of the interface between soil layers were derived from three types of images, namely the true electrical resistivity model, an ERT tomography, and an EMI tomography. The geometry of the interface was estimated for each type of image with the potential field method. Finally, the estimated geometry was compared to the true geometry of the synthetic model.
Each layer of the synthetic model was simulated as a realisation of a Gaussian random field via covariance matrix decomposition (Constantine, 2020). The random fields were characterised by an anisotropic covariance model where t x and t z are the horizontal and vertical distances between the pair of points x 1 and x 2 , s 2 is the variance of the random field, and u x and u z are the horizontal and vertical correlation lengths, respectively. The parameters of the random fields are summarised in Table  1. The air-ground interface was considered flat and located at z = 0 m. The interface between the soil layers was located at In this case, the position of the interface also corresponded to the thickness of the upper layer. The synthetic model served as the base for simulating geophysical data and also as a representation of a tomography with perfect resolution. ERT and EMI tomography were constructed by simulating data acquisition on the synthetic model. The ERT tomography was constructed with data simulated following a roll-along pattern. The EMI tomography was constructed with EMI data simulated with the acquisition geometry of the commercial device DUALEM-421 (Section A.2). The simulation parameters for the ERT and EMI data sets are summarised in Table 2.  The Laplacian and gradient magnitude sufficed to detect the interface between soil layers in the true model, so smoothing was not applied. Figure 3(b) shows the edge detection technique applied to the original ERT tomography. Visually, the ERT tomography indicates the presence of two soil layers with internal variability. The zero-crossings of the Laplacian showed both true soil interfaces and fake edges. The fake edges arised from small-scale variability in the ERT tomography. Smoothing and thresholding were applied to improve the detection of edges in the tomography. Figure 3(c) shows the ERT tomography after a Gaussian filter with a standard deviation of 2.5 was applied. Although the filtered ERT tomography was not visually different from the original tomography, the small-scale zerocrossings were significantly fewer. Two geometrically similar edges were visible in the Laplacian of the filtered tomography. The edge that corresponds to the soil interface was identified by thresholding the gradient magnitude. Figure 3(d) shows the edge detection technique applied to the EMI tomography. Visually, the EMI tomography also indicates the presence of two soil layers. However, the contour plot of the unfiltered EMI tomography resembles only roughly to the true interface geometry. Moreover, the Laplacian of the unfiltered tomography is contaminated with small-scale zerocrossings, so the interface between soil layers is not visible. Thus, smoothing was applied. Figure 3(e) shows the EMI tomography after a Gaussian filter with a standard deviation of 20.0 was applied. The remaining edge in the Laplacian is located at the same position where the gradient magnitude is large.
3.3. Potential field: true model, ERT, and EMI Figure 4(a) shows the potential field method applied to the data sets of the true model, the ERT, and the EMI tomography. The data sets consist of contact points and orientations. The contact points between the soil layers were obtained from local point data of the subsurface, i.e. Equation (8) sampled at x=23 m and x=46 m. The orientations were calculated as the dip angle of the Laplacian interfaces. The Laplacian interfaces are the edges detected with the Laplacian edge detection technique (Figure 3). They represent the interface between soil layers estimated with geophysical data only. The potential-field interfaces represent the interfaces estimated by combining geophysical data and local point data of the subsurface. The potential-field interfaces correspond to one of the contour lines of the potential field. Figure 4(b) shows the error incurred by the Laplacian and the potential-field interfaces. The errors of the Laplacian and potential-field interface of the true model are negligible. Thus, sampling the subsurface is not necessary in the hypothetical case that the image of the subsurface is perfect. The added value of a combined approach for site investigation is visible when the image of the subsurface is approximate. In the ERT and EMI data sets, the error of the potential-field interface was smaller than that of the Laplacian interface. Thus, the combination of geophysics and local point data improved the estimation of geometric variability.

Study site and data collection
The case study site is located in the central-northern section of the Rhine-Meuse delta. In this alluvial environment, multiple paleochannels are present in the subsurface. These channels are visible in the present-day landscape as topographic ridges ( Figure 5). Within this area, the Stuivenberg channel belt is the main topographic expression. Berendsen (1982) investigated comprehensively the geological development of the area. Additionally, Winkels et al. (2017) used local borehole data to gain insights into the internal buildup of the Stuivenberg channel belt and encasing sediments. Figure 5 shows an approximate lithological cross-section from Winkels et al. (2017). The cross-section shows the Stuivenberg paleochannel which mainly consists of sandy deposits. The channel is surrounded by clayey deposits and a Pleistocene sandy substrate. The geometric variability of the clay-sand interface was investigated at the location of the lithological crosssection. Since the clay layer is in contact with air, the depth of the clay-sand interface also represents the thickness of the clay layer.
The clay-sand interface was identified in boreholes from Winkels et al. (2017). The data set consisted of eight borehole cores (Table 3), each of which was sampled with a resolution of 0.1 m. The interface depths were determined by visual inspection. Tomographic images, collocated with the lithological cross-section, were created with ERT and EMI data. The ERT data were collected in spring 2018 with a Wenner-alpha array. The roll-along technique was used to cover the extent of the test site. The EMI data were collected in summer 2020 with a DUALEM-421 (Section A.2). Table 4 summarises the acquisition parameters for the ERT and EMI data. Clayey and sandy layers were recognised in ERT and EMI tomography by their value of electrical resistivity. In general, low values of electrical resistivity correspond to clayey layers while high values correspond to sandy layers.

Edge detection and interpretation: ERT and EMI tomography
Figure 6(a) shows the Laplacian edge detection technique applied to the ERT tomography. The tomography was filtered with a Gaussian filter at a standard deviation of 4.0. The value was chosen so that near-surface artefacts were reduced without distorting the tomography. Due to the variability of the geological setting, multiple soil layers were detected in the Laplacian of the ERT tomography. The interfaces that correspond to the layer of interest, the bottom of the clay layer in Figure  5, were picked from the edges detected in the tomography. The picked edges showed a large value of gradient  magnitude. Figure 6(b) shows the Laplacian edge detection technique applied to the EMI tomography. The tomography was filtered with a Gaussian filter at a standard deviation of 10.0. The value was chosen so that near-surface artefacts were reduced without distorting the tomography. Although the EMI tomography does not resolve the geological setting in the same detail as the ERT tomography, similar insights were derived.
The clayey deposits on top of the Stuivenberg sand body (from approximately 60 to 350 m in Figure 5) are thinner than the deposits on the sides of the sand body. On top of the sand body, the thickness of the clay layer varies at a smaller scale due to autogenic processes. For instance, BH1 and BH2 (Table 3), which are located on top of the sand body, show a thicker clayey layer than the surrounding boreholes. Small-scale variations are visible in both ERT and EMI tomography. However, only the Laplacian of the ERT tomography captures the details of the fine-scale variability. The Laplacian of the EMI tomography does not capture this variability due to the lack of detail of the tomography.
The lack of detail in the EMI tomography is inherent to the device used in this survey, DUALEM-421. Two factors affected the EMI tomography. First, the measurements were collected during a dry summer day. Therefore, the electrical resistivity of the shallow portion of the ground was high due to evapotranspiration. Higher electrical resistivity values resulted in lower signal strength and thus lower resolution. Second, the maximum depth of investigation of the DUALEM-421 is approximately 6 m. Meanwhile, the thickness of the clay layer is known to be larger than 6 m on the left-hand side of the sand body. Therefore, the claysand interface in the EMI tomography was only roughly visible on this side. On the right-hand side of the sand body, the thickness of the clay layer is known to be smaller than 4 m. Therefore, the clay-sand interface is more visible on this side.

Potential field: ERT and EMI
Figure 6(g) shows the potential field method applied the ERT data set of the Montfoort site. The estimated geometry of the clay-sand interface corresponds to one of the contour lines of the potential field. The method was applied with borehole data and ERT tomography. The location of the clay-sand interface was derived from borehole data (Figure 6(c)) and the orientations were derived from ERT tomography (Figure 6(e)). The estimated interface with the ERT data set captures in great detail the complexity of the lithological cross-section of the Stuivenberg channel belt.
The geometry of the clay-sand interface was also estimated with the EMI data set. For that purpose, additional information was manually added regarding the orientation of the clay-sand interface on top of the sand body. Although finer variability is visible in the EMI tomography, the Laplacian edge detection technique was not successful at detecting variability in this region. Thus, we considered the clay layer on top of the sand body to be horizontal (Figure 6(f)). The orientations were positioned at the minimum depth found in the boreholes from Table 3 which is 0.6 m. Additionally, the data from a borehole core which is not listed in Table 3 was used (BH9 in Figure 6(d)). In the borehole core, which was drilled at x=30 m, he Pleistocene sand layer was not reached after 5.5 m of drilling. Thus, the location of the clay-sand interface was conservatively estimated at a depth of 5.5 m. Figure 6(h) shows the potential field method applied to the EMI data set. To a lesser detail than the ERT tomography, the clay-sand interface was estimated with the EMI data set. On average, the clay-sand interface estimated with the ERT data set is comparable to that with the EMI data set (Table 5).

Study site and data collection
The site is a dike stretch which is part of the dike trajectory 20-3. The trajectory 20-3, which goes along with the Spui river, is part of the dike ring 20 Voorne-Putten in the Netherlands. The site is located in a tidal environment where tidal paleochannels are expected, but not  Figure 7 shows the elevation map of the study site and an approximate lithological cross-section. The cross-section shows a roughly homogeneous upper clayey layer underlain by a sandy Holocene layer. An isolated tidal paleochannel crosses both layers. The geometric variability of the clay-sand interface was investigated at the location of the lithological cross-section. As in Section 4, the depth of the clay-sand interface also corresponds to the thickness of the clay layer.
Local authorities have carried out extensive CPT investigation along the dike trajectory 20-3. Two CPTs were available at the location of the study site. The depth of the clay-sand interface was obtained from the Soil Behavior Type (SBT) index of the CPTs (Robertson, 2009). The interface between the upper and lower layer was defined the depth were the SBT index reached a  value of 2.6 (Robertson, 2009). Tomographic images, coincident with the lithological cross-section, were created with ERT and EMI data. The ERT data were collected with a Wenner-alpha array. The roll-along technique was used to cover the extent of the test site.
The EMI data were collected with a DUALEM-421 (Section A.2). Table 6 summarises the acquisition parameters for the ERT and EMI data. Clayey and sandy layers were recognised in ERT and EMI tomography by their value of electrical resistivity. In general, low values of electrical resistivity correspond to clayey layers while high values correspond to sandy layers.

Edge detection and interpretation: ERT, and EMI tomography
Figure 8(a) shows the Laplacian edge detection applied to the ERT tomography of the study site. The tomography was not filtered because small-scale artefacts were not strongly present. Multiple layers were detected in the Laplacian of the ERT tomography. The tomography shows a high resistivity top which is the dry part of the upper layer. After the dry part, the upper layer shows low resistivity. The layer after shows higher resistivity values. The variability of the clay-sand interface is visible in the Laplacian of the ERT tomography (Figure 8 (a)). Figure 8(b) shows the Laplacian edge detection applied to the EMI tomography of the study site. To reduce small-scale artefacts, the tomography was filtered with a Gaussian filter at a standard deviation of 10.5. Although the depth of the EMI tomography is shallower than that of the ERT tomography, both tomographic images show comparable features. The tidal paleochannel of Figure 7, which is located from x=50 m to x=100 m, was captured in the ERT and EMI tomography. In the ERT tomography the paleochannel was fully captured while in the EMI tomography only the top was captured. The clay-sand interface is visible in the Laplacian of the ERT tomography along the entire dike stretch. Meanwhile, the Laplacian of the EMI tomography shows the clay-sand interface mostly at the location of the paleochannel intrusion. The penetration depth of the EMI survey device, DUALEM-421, was lower than the depth of the clay-sand interface on the sides of the paleochannel.

Potential field method with ERT and EMI
The potential field method was applied to the ERT and EMI data set. The CPTs were located at the top of the sand body intrusion and at the flat portion of clayey layer (Figures 8(c,d)). The clay-sand interface in CPT0 was located at a distance of 78.0 m and a depth of −1.5 m while the clay-sand interface at CPT1 was located at a distance of 178.0 m and a depth of −6.0 m. A fictitious interface at a distance of 25.0 m and depth of −6.0 m was added to the EMI data set (CPT2 in Figure 8(d)). The interface was added to enforce the clay-sand interface at this location to be similar to the interface at 178.0 m. The orientations derived from the ERT and EMI tomography are shown in Figures 8 (e) and 8(f). Since the depth of the EMI tomography captures the top of the sand body intrusion, but it does not reach the deeper sand layer. We manually added the orientations of the layer interface based on the interpretation of neighbouring CPTs. We considered this layer to be horizontal. The orientations  derived from the EMI tomography and the manuallyadded orientations are shown in Figure 8(f) in black and blue respectively. The estimated geometry of the upper clayey layer is shown in Figure 8(g) for the ERT data set and in Figure 8(h) for the EMI data set. The geometric variability of the clay-sand interface is captured by both data sets. On average, the clay-sand interface estimated with the ERT data set is comparable to that estimated with the EMI data set (Table 7).

Discussion
In the Netherlands, stochastic subsurface models were developed for the reliability assessment of primary  dikes (Hijma and Lam, 2015). The models were created based on an extensive database of local point data and geological knowledge of the Rhine-Meuse delta. Despite the large amount of available information in the Netherlands, uncertainties remain in the subsurface models. The main source of uncertainty is the sparsity of the data relative to the geological variability. To account for uncertainty in reliability assessments, geological experience is often used. For instance, a probability of occurrence is assigned to a geological scenario even though that scenario is not found in site investigation data (ENW, 2012). Such geological scenarios could lead to over conservative reliability assessments. To discard or confirm this type of geological scenarios, it is necessary to reduce data sparsity. Geophysical exploration methods are a powerful tool in this regard because they provide a horizontally-continuous view of the subsurface. In this research, the geometric variability of soil layers was studied. The geometry of soil layers is an important part in the reliability assessment of dikes especially for the failure mechanisms of piping and macrostability. The approach presented here aims at improving the characterisation of geometric variability in two ways. First, by showing that geophysical exploration can improve the allocation of site investigation efforts. Second, by improving the estimated geometry of soil layers with a reduced amount of local point data.
Geophysical exploration can improve the allocation of site investigation efforts by providing a view of the average composition of the subsurface. For example, the geophysical images in Sections 4 and 5 (Figures 6  and 8) show the regions of the subsurface where the geometry of soil layers varies due to the presence of paleochannels. Indeed, anomalies, such as paleochannels, are a large source of uncertainty in reliability assessments because they are easily missed in local point data (ENW, 2012). Elevation maps are a valuable source of information to visually detect paleochannels (Berendsen and Volleberg, 2007). However, small paleochannels, e.g. Section 5, are more challenging to detect in elevation maps. Then, geophysical methods can fill the information gap between local point data. A common problem of geophysical exploration is that the data interpretation is often subjective and relies on expert knowledge. Thus, an automatic edge detector was applied to the geophysical images so that the interpretation is objective and reproducible without expert knowledge. The information retrieved from geophysical data is valuable for assisting the allocation of site investigation efforts. For example, by indicating the location and extent of anomalous features, such as paleochannels.
Apart from assisting exploration efforts, geophysical data improves subsurface characterisation when combined with local point data. In particular, the estimation of geometric variability is improved. To elaborate, we consider the estimation error of the ERT data set of the Montfoort site. Figure 9 shows the estimation error of the Montfoort site. The left-hand side of the plot shows the absolute error per borehole. The righthand side shows the overall density distribution of the Montfoort site. The absolute error in Figure 9 was calculated using the contact points measured in boreholes (Table 3) as the reference truth. Figure 9(a) shows the error incurred when only the Laplacian interface of the ERT tomography is used for the interface estimation. In this test site, there is a good agreement between the Laplacian and borehole interface. Indeed, the density distribution of the error is narrow with few outliers . Figures 9(b) and 9(c) show that the estimation of geometric variability improves when geophysical and borehole data are combined via the potential field method. The error in Figure 9(b) was calculated by removing one borehole at a time from the estimation procedure and using that borehole for validation, namely cross-validation. The median of the error in Figure 9(b) is reduced significantly and the overall density distribution of the error is narrow. The estimation with a combined data set shows consistency even when the number of boreholes is reduced. Figure  9(c) shows the cross-validation error when half of the boreholes are used in the estimation. The figure shows that the median of the error remains low and the density distribution remains narrow. Figure 9(c) shows more data points because there are several possible combinations to remove four boreholes.

Conclusions
An approach to estimate the geometry of soil layers with the aid of geophysical tomography was presented. The approach aims to address the sparsity problem existing in the reliability assessment of dikes, i.e, the mismatch between site investigation density and geological variability. The approach provides a two-fold quantitative framework to incorporate geophysical information into geological schematisation. First, an automatic edge detection technique allows for objective and reproducible interpretation of geophysical exploration data. Second, a cokriging framework allows for objective incorporation of geophysical and local point data to estimate the geometric variability of soil layers. The resulting estimation is improved even when a reduced amount of local point data is available. In addition to this, geophysical data showed to be a valuable tool for assisting the allocation of site investigation efforts.
The approach presented in this study compared the data from two geophysical methods, namely ERT and EMI. The ERT method showed a detailed description of the subsurface. To a lesser detail than the ERT method, the EMI method also showed the main features of the subsurface, such as the presence of paleochannels. From an operational point of view; however, the EMI method is a more viable method for exploring the subsurface of dikes, for it can cover wide areas in a fraction of the time required with the ERT method.
Looking at the geophysical data beyond the geometry of soil layers, it is clear that a larger degree of complexity is captured in the geophysical data. The challenge is to retrieve quantitative information from this data for geotechnical applications. Geophysical insights are needed in terms of, for example, the internal correlation structure of soil layers.
De La Varga from Aachen University provided us with insights into the numerical aspects of the potential field method.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
This work is part of the research programme All-Risk with project number PP15-21, which is (partly) financed by the Dutch Research Council (NWO).

A.1. Electrical resistivity tomography forward model
The ERT method measures the ground response, electric potential, of point sources of direct current. The injection and the measurement electrodes are commonly located at the surface of the medium. However, they can also be located inside the medium. In the case of stainless steel electrodes, the injection and measurement points are interchangeable. Thus, an electrode can be used to inject current or to measure electric potential. The ground response depends on the bulk electrical resistivity of the ground. For a point current source, the ground response is related to the bulk electrical resistivity by the equation where r(x, y, z) is the electrical resistivity of the subsurface, f is the electric potential, I is the source current, and δ is the Dirac-delta function centred at the location of the point source. Equation (A1) represents the forward model for ERT. In reality, two sources of opposite polarity are needed to inject current into the ground. Also, two electrodes are needed to measure the electric potential. Thus, one ERT measurement requires four electrodes, a quadrupole. Figure  A1 shows a schematic of the ground response in an ERT measurement. In the figure, C 1 and C 2 are current electrodes while P 1 and P 2 are potential electrodes. The total electric potential results from the algebraic sum of the electric potential for each point source. An ERT measurement is represented by the apparent resistivity instead of the electric potential. The apparent resistivity is the resistivity of an equivalent homogeneous medium. In a homogeneous half-space, the ground response for a single point source is given by where r is the distance between the current and potential electrode. Thus, the electric potential difference between P 1 and P 2 becomes where r C i P j is the distance between injection point C i and measuring point P j . Thus, the apparent resistivity becomes where To detect different portions of the subsurface, the separation between electrodes or the array type is changed. Commonly used array types are dipole-dipole, Wenner-alpha, and Wenner-Schlumberger. The array type has a large influence on the sensitivity of the measurement. For example, the Wenneralpha array is more sensitive to horizontal structures while the dipole-dipole array is more sensitive to vertical structures (Loke, 2013). The python package pyBERT (Rücker, Günther, and Spitzer, 2006) was used to solve the forward model in equation (A1).

A.2. Electromagnetic induction forward model
The EMI method described in this section measures the ground response, magnetic field, of a vertical magnetic dipole above the earth surface. The magnetic dipole is induced by a coil through which alternating current flows at a given frequency. The geometric arrangement of the source and receiver coil with respect to each other determines the sensitivity of the signal. The commercial device DUALEM-421 (Dualem Inc., Milton, ON, Canada) was used in this study. The geometric arrangement of this device consists of three horizontal coplanar (HCP) sensors and three perpendicular coplanar (PRP) sensors. Thus, each measurement contains six data points. The separation between transmitter and receivers is fixed. The HCP sensors are separated from the transmitter by 1.0, 2.0, and 4.0 m. The PRP sensors are separated from the transmitter by 1.1, 2.1, and 4.1 m. Figure A1. Ground response, electric potential f, of a homogeneous medium under two point sources of opposite polarity. Current I = 1 (A) and electrical resistivity r = 1 (ohm m). C 1 and C 2 are current electrodes. P 1 and P 2 are potential electrodes. Figure A2 shows the ground response of the vertical magnetic dipole, m, induced by a DUALEM-421. The magnetic field measured at the receiver coils consists of the primary and secondary magnetic field. The primary field is that measured in free space while the secondary field is that reflected by the ground. For a vertical magnetic dipole, the magnetic field has a radial component (H r1m , H r2m , and H r4m in Figure A2) and a vertical component (H z1m , H z2m , and H z4m in Figure A2). The vertical component of the magnetic field is measured by the HCP sensors while the radial component is measured by the PRP sensors. At the receiver coils, the secondary magnetic field is calculated only from the upgoing component of the potential F − 0 e −u0h e u0z . A detailed derivation of the solution of Maxwell's equations for a vertical magnetic dipole in a layered medium is presented in Ward and Hohmann (1988). In this section, we present only the solutions. The primary magnetic field is given by where m is the magnetic dipole strength and s is the distance between transmitter and receiver. The measured response given by commercial devices is given as the ratio of the secondary magnetic field to the primary magnetic field H s /H p . Thus, the strength of the magnetic dipole cancels out. The ratio H s /H p for the HCP sensors is given by Meanwhile, the ratio H s /H p for the PRP sensors is given by In the above equations, r TE is the reflection coefficient, λ is the wavenumber, and J i is the ith-order Bessel function of the first kind. The reflection coefficient, is a function of the wavenumber,û 1 , which is calculated iterativelyû 1 = u1û 2 + u 1 tanh (u 1 h 1 ) u 1 +û 2 tanh (u 1 h 1 ) (A10) u n = unû n+1 + u n tanh (u n h n ) u n +û n tanh (u n h n ) where m n and s n are the magnetic permeability and electric conductivity of layer n and w is the angular frequency of the transmitted signal. In many applications, the magnetic permeability is equal to that of free-space. Therefore, the quantity that is mapped by EMI is the electrical conductivity of the ground. Equations (A7) and (A8) are the forward model for EMI of a DUALEM-421.

A.3. Tomographic inversion
A geophysical tomography is constructed from geophysical data via inversion. The tomography represents the spatial distribution of material properties in the subsurface. Inversion consists in minimising the squared error between the measured data and the simulated data from a forward model. The error is minimised by changing the material properties of the medium. In the case of ERT, the forward model is Equation (A1) from which the simulated data are the apparent resistivity, Equation (A4). For EMI of a DUALEM-421, the simulated data are given by Equations (A7) and (A8). The inverse problem for ERT and EMI is ill-posed. Thus, mathematical constraints, regularisation, are needed to find a unique solution. A tomography is constructed by minimising an objective function, Φ, which contains a data error component and a regularisation component The first term in the above equation contains the data error component in which d is the measured data, f (m) is the simulated data with the forward model, f, applied to the material properties of the medium m, D is a weighting matrix that contains the correlation structure of the measurement error. In the second term λ is the regularisation strength, m 0 is an a priori model, and C is the constraint matrix that contains information about the expected correlation structure of the final model. Smoothness-type of constraints are commonly used in geophysical inversion. The tomography, m, is the model that minimises the objective function from equation (A15). A comprehensive treaty of inverse theory is presented in Menke (2012). The python package pyGIMLi (Rücker, Günther, and Wagner, 2017) was used for tomographic inversion.

A.4. Potential field method
The covariance model and drift function model, which are needed to formulate the cokriging system of Equation (5), are elaborated in this section. To aid the illustration, a typical data set consisting of contact points and orientations is used ( Figure A3(a)). In Figure A3(a) and the following equations, the contact point i associated to layer p is denoted as x pai Figure A3. Potential field method. (a) Example data set used for the application of the potential field method. x pai indicates the contact point i associated to layer p. x bk indicates the orientation point k. (b) Potential field (contour lines) estimated from the data in (a). The contact points from the same layer share the same contour line.
while the orientation point k is denoted as x bk . The available data consist of the random function Z 1new evaluated at the contact points x pai and the random function Z 2 evaluated at the orientation point x bk , i.e, the orientation at point x bk . The covariance of Z 1new is defined as an expectation which simplifies to where x pa0 and x qa0 are the reference points for layer p and q, respectively. The covariance of Z 1 is defined in terms of the absolute distance, r, between data points. Meanwhile, the covariance of Z 2 is defined in terms of the gradient components Z 2x , Z 2y , and Z 2z . Following Lajaunie, Courrioux, and Manuel (1997), the covariance between the gradient points x bk and x bl becomes and where ru and rv are the signed distances between points along the components x, y, or z. The cross-covariance between the contact point x pai and the gradient point x bk becomes The drift functions of Z 1new and Z 2 describe the trend of the random functions. The drift function of Z 1 can be described with a polynomial. For example, a polynomial of first degree U Z1 (x pai ) = 1 + x + y + z (A21) Figure A4. Generalised apparent resisitivity response for a two-layer model and a Wenner-alpha array. Figure A5. Particular apparent resisitivity response for a two-layer model and a Wenner-alpha array.
or a polynomial of second degree U Z 1 (x pai ) = 1 + x + y + z + x 2 + y 2 + z 2 + xy + xz The drift of Z 1new follows from the drift of Z 1 The drift of Z 2 is calculated as the gradient of U Z 1 where u is the component x, y, or z. The covariance and cross-covariance matrices require a covariance model for Z 1 . However, the covariance of the potential field cannot be derived experimentally because the potential field is not known. On the other hand, the experimental covariance of the gradient Z 2 can be derived from the data because the orientations are known (Aug, 2004). Alternatively, a heuristic covariance model can be used.
and a default value for the range where L, B, and H are the dimensions of the model domain. A detailed justification for such model and values is presented in Chiles et al. (2004). Figure A3(b) shows the estimated geometry using a cubic covariance model and the default values for the variance, C 0 , and range, a. The plot shows that the contact points of the same layer belong to the same iso-surface of the potential field. The interpolation of Figure A3(b) was generated with a python implementation of the potential field Figure A6. Edge detection applied to one-dimensional tomographic images before and after a smoothing Gaussian filter was applied.
method. The implementation is added as supplementary material.
Appendix 2. Edge detection to detect soil layers: Two-layer model A synthetic two-layer model is studied with the ERT method. First, the impact of geophysical measurements in a tomographic image is studied. For that purpose, the analytical ERT response for a generic two-layer model is presented. Subsequently, the study is extended to particular cases of soil properties and acquisition designs. The synthetic responses from these particular cases are used to construct one-dimensional tomographic images. Finally, the performance of the Laplacian edge detection on these images is studied. Figure A4 shows the apparent resistivity response for a Wenner-alpha array in a horizontal two-layer model. The apparent resistivity is normalised with respect to the resistivity of the upper layer. The apparent resistivity contour is plotted with respect to the normalised electrode spacing, a/h, and the resistivity contrast k r = r 2 − r 1 r 2 + r 1 .
The normalised electrode spacing, a/h, determines the sensitivity of the measurement in depth. A measurement with a small value of a/h is more sensitive to the resistivity of the upper layer. Meanwhile, a measurement with a large value of a/h is more sensitive to the lower layer. An acquisition array should contain small and large values of a/h so that the upper layer and lower layers are properly detected. The optimum range of a/h values will depend on the resistivity contrast, k r . At a given contrast (white dashed lines in Figure  A4), the range of a/h values has to be such that the resistivity response is properly sampled. In other words, an acquisition array should cross the contour lines in Figure A4 from r a /r 1 = 1 to r a /r 1 = r 2 /r 1 . Operational limitations regarding data acquisition and lack of geological knowledge, make comprehensive data acquisition difficult to achieve in practice. Thus, the measured data contain incomplete geological information which affects the tomography. A particular two-layer model with h=1.5 m and r 1 = 10 ohm m is further analysed. Two resistivity contrasts are analysed, namely k r = 0.2 and k r = 0.6 which result in r 2 = 15 ohm m and r 2 = 40 ohm m, respectively. Figure A5 shows the noise-free apparent resistivity response for each resistivity contrast. For tomographic inversion, four cases are considered each with 23 measurements. The base electrode spacing, a, for each case is 0.25, 0.5, 1.0, and 3.0 m. Due to the limited number of measurements, the apparent resistivity response is not fully reconstructed in the measurements. The portion of the response that is covered in each case is shown at the bottom of Figure A5. At a = 0.25 m, the response is more sensitive to the upper layer than the lower layer. At a = 0.5 m and a = 1.0 m, the apparent resistivity is sensitive to both layers. Even though the measurements do not reconstruct the complete apparent resistivity response, the trend of the response is well captured. At a = 3.0 m, most of the measurements are sensitive to the lower layer while few measurements are sensitive to the upper layer.
The synthetic apparent resistivity data were contaminated with Gaussian random noise of 2% before tomographic inversion. The value follows the guideline of the report for Best Practices in Electrical Resistivity Imaging (Day-Lewis et al., 2008). Figure A6 shows the Laplacian edge detection technique applied to the tomographic images obtained for k r = 0.2 and k r = 0.6. The figure shows the tomographic images and edge detection before and after a Gaussian filter was applied. The chosen standard deviation for the filter was 5.5 ohm m. Most of the inverted models reproduce the trend of variability of the true model. However, the inverted models at low resistivity contrast, k r = 0.2, are more affected by the added noise than the inverted models at high resistivity contrast, k r = 0.6. This fact is reflected in the undulations of the inverted models which result from overfitting the data during tomographic inversion. In general, a combination of low resistivity contrast and nonoptimum data acquisition is detrimental to the reliability of tomography for geological interpretation. However, no exact limit can be drawn. For example, the inverted model for k r = 0.6 and a=3.0 m shows a clear mismatch between the true and inverted resistivity values even though the contrast is high. On the other hand, all the inverted models for k r = 0.2 reproduce the trend of electrical resistivity. A sensible criterion to assess the reliability of a tomographic image is that layers that are visible in the tomography should also be visible in the raw data. In other words, the measured response should match that of a conceptual lithological cross-section. Figure A6 shows Laplacian edge detection technique applied to one-dimensional tomographic images for k r = 0.2 and k r = 0.6. The automated edge detector finds edges in the tomographic images even when there are no geological interfaces. The so-called fake edges are more often found when smallscale artefacts are present in the tomography. For instance, fake edges are more frequent for k r = 0.2 because of the undulated pattern of the inverted models. In Figure A6, a smoothing Gaussian filter reduced the undulations of the tomographic images significantly and therefore fake edges. At the same time, the filter kept intact the large-scale trend of resistivity values. Fake edges were further reduced by setting a threshold on the gradient. The tomography and edges detected after filtering, and thresholding agree with the resistivity trend and the interface location of the true model. However, the locations of the tomographic interfaces do not precisely coincide with that of the true model. Hence, the potential field method (Lajaunie, Courrioux, and Manuel, 1997) was used to calibrate the tomographic data with more accurate information obtained from direct samples of the subsurface.