## Abstract

Stratigraphic sequences in boreholes are commonly estimated by interpreting combinations of well logs. The interpretation is generally tedious and is made some time after log completion, which often leads to a loss of valuable first-hand information gathered on-site. This may lead to delayed or potentially poor on-site decisions. To make things worse, the standard interpretation of well logs is, at least to a certain degree, subjective and based on the manipulation of data, which may be difficult to trace in the long term. Small changes in lithology are often disregarded and alternating thin layers presenting different lithologies are often combined in one single (notably thicker) stratigraphic unit. Therefore an automatic parameter-based and thus traceable and objective quick look at the lithology immediately after log completion represents both a valuable tool to help with on-site decisions and a solid, mathematically based starting point for further physically based interpretations carried out by log analysts. We present a workflow for the interpretation of well logs defined as an optimization problem. The workflow is applied to the characterization of metre- to decametre-scale stratigraphic units along 13 boreholes in northern Switzerland (one-dimensional resolution) and to millimetre-scale features over a wall at the Mont Terri underground rock laboratory in Switzerland (two-dimensional resolution). The results show that: (1) the workflow accurately maps lithological changes; (2) the interaction with the analyst is minimized, which reduces the subjectivity of the interpretation; and (3) outputs are available for on-site decisions.

Geophysical well logs are one of the best sources for identifying facies interfaces or characteristic hydrogeological properties in boreholes (Dewan 1983; Crain 1986). Thus it is not strange that well logging has become a standard in the oil, gas and geothermal industries and in hydrogeology. Facies mapping, the identification of formation fluids, correlation between boreholes and evaluation of the productive capabilities of reservoirs are the main objectives of the interpretation of well logs (Murray & Geldart 1991).

Facies mapping (in a broad sense, the estimation of a lithostratigraphic sequence) is generally carried out by the visual inspection of a composite of well logs or suitable transformations of them (Crain 1986; Rider 1986). As such, standard interpretations are based on subjective criteria (e.g. filters or thresholds), which may lead to different analysts recognizing different lithostratigraphic sequences. In addition, the transformations required for such analyses are difficult to trace in the long term. In addition to subjectivity, the standard interpretation of well logs presents two additional drawbacks. First, small changes in lithology (e.g. in the clay content) are often disregarded and alternating thin layers presenting slightly different lithologies are often combined in one single, notably thicker, lithostratigraphic unit. This is not a problem if the goal of the interpretation is to capture large-scale trends. However, small-scale features may be clearly observable in geophysical well logs (the sampling interval of which ranges from 0.5 to 15 cm according to Schlumberger standards) and their accurate identification and mapping may be crucial because they may be preferential flow paths for contaminant migration or water-conducting features. Second, the interpretation is generally made some time after log completion, often leading to a loss of valuable first-hand information gathered on-site (e.g. during drilling). Therefore an automatic, parameter-based and hence objective and traceable quick look at the lithology immediately after log completion represents a valuable tool that helps with on-site decisions (e.g. the depth of hydraulic tests).

A lot has been achieved since the first resistivity log was recorded in 1927 by Schlumberger (at Pechelbronn field in Alsace, France; Clapp 1932): (1) to enhance the logging accuracy to capture small-scale features (although these are generally lost during standard log interpretation); (2) to reduce the subjectivity in the interpretation process; and (3) to facilitate the interpretation by increasing its speed, making it beneficial for on-site decisions. The first algorithm to automatically identify facies boundaries is attributed to Lanning & Johnson (1983) and this has been successfully applied numerous times (Ben-Hamadou *et al.* 2005; Maiti & Tiwari 2005; Mokhtar 2007; Al-Garni 2008; Xiang *et al.* 2009; Javid & Tokhmechi 2012; Ofuyah *et al.* 2014; Singh *et al.* 2016). The algorithm is based on a Walsh transform of the recorded well logs. A Walsh transform is analogous to the Fourier transform with a constant window for analysing stationary signals (Pan *et al.* 2008). Window size is addressed by the short-time Fourier transform (Allen 1977; Allen & Rabiner 1977). The main drawback is that only one window (narrow for high-frequency data or wide for low-frequency data) can be used at a time. The windowing problem can be tackled by wavelet transforms that accommodate different windows for different frequency bands (Nason & von Sachs 1999; Polikar 1999). The set of windows, and therefore the discrete wavelet transform, is estimated by thresholding techniques. Different thresholds can be applied (leading to different results; Nason 1995), which does not help to decrease the subjectivity of the interpretation. In addition, wavelet transform methods are not easy to implement and are not as intuitive as the original method suggested by Lanning & Johnson (1983).

Regardless of the transformation applied to enhance transitions in well logs (e.g. Walsh, short-time Fourier or wavelet transform), the mapping of lithological interfaces requires the definition of a set of parameters that define the mapping accuracy and a set of weights balancing the quality and quantity of the information contained in each log. The choice of such parameters has an impact on mapping and is subjective, tedious (usually by trial and error) and difficult to trace in the long term. This work is aimed at filling these gaps.

We present a new workflow for the automatic interpretation of well logs based on the Lanning and Johnson algorithm that: (1) is parameter-based and therefore traceable; (2) is fast, enabling on-site decisions based on well logs and other sources of information; (3) is accurate because it is based on actual data and therefore can resolve small-scale features if these are captured by the well logs; and (4) accommodates the calibration of the intervening parameters, substantially reducing the subjectivity of the interpretation process. Calibration requires some prior information of what is being modelled, e.g. the depth of facies interfaces. Such information may arise from two non-exclusive sources. It can be provided directly on-site (e.g. by the driller who perceives small fluctuations in the rate of penetration or by the field geologist who monitors the cuttings brought to the surface) and/or, if a posterior analysis is made, information from core inspections/sampling or; the interpretation made by a log analyst may serve as calibration data. Both are valuable. On-site information can be as accurate as the analyst's interpretation and often yields valuable information on the small-scale variability, which is generally lost after standard interpretation of well logs. Even if a first (subjective) interpretation has been made, it is convenient to recall the original datasets in an objective manner because important small-scale features disregarded in the first analysis may be revealed.

This paper is organized as follows. First, the workflow is outlined and illustrated using well logs in the borehole BDB-1 at the Mont Terri Underground Rock Laboratory (Mont Terri URL, Switzerland). Second, the workflow is applied to mapping large-scale lithostratigraphic units in 13 boreholes in northern Switzerland and to the characterization of two-dimensional millimetre- to centimetre small-scale features at a wall in the Mont Terri URL. The paper finishes with some conclusions and recommendations about the use of automatic tools for the objective, traceable and quick interpretation of well logs.

## Workflow

The workflow consists of four main blocks (Fig. 1). In block 1, automatic transformations of well logs are carried out to make them comparable. In block 2, rock boundaries (in a broad sense, lithofacies interfaces) are automatically detected. In block 3, the statistics of the well logs or correlated variables (e.g. shale volume, porosity or hydraulic conductivity) are calculated at each recognized lithostratigraphic unit. Note that blocks 2 and 3 are embedded in the calibration loop. The last block of post-processing allows the calculation of variograms, which enables further geostatistical modelling – for example, the correlation of logs at different boreholes. The different steps are illustrated in the following sections of the paper by applying the workflow and comparing the results with values taken from a lithostratigraphic profile of borehole BDB-1 at the Mont Terri URL (stratigraphic profile from Reisdorf *et al.* 2014; geophysical data obtained from the Mont Terri Consortium).

### Block 1: prior transformations

The workflow allows the consideration of several well logs simultaneously. These are sometimes collected with different logging tools and measurement frequencies or sampling intervals (Fig. 2). In addition, the scales of the logs and even the physical units may be different (e.g. density log may be presented in kg m^{−3}, but the gamma logs in API units), making them hardly comparable. To alleviate this problem, the first step consists of introducing a common discretization and scale for all logs. Each log is resampled at the points defined by the chosen discretization (usually that of the log with best resolution to avoid loss of information) by means of a linear interpolation scheme and is standardized to its maximum value (therefore all in the range [0, 1]). The discretization is a filtering parameter, which needs to be calibrated.

Borehole quality also plays an important part in the mapping because some logs (e.g. gamma ray logs) are not representative near or at borehole breakouts because the rock/formation volume under analysis is smaller than that in the absence of breakouts. Such borehole singularities are easily detected by observing the caliper log. Measured signals negatively impacted by borehole quality are reconstructed along these intervals by a multi-point statistical algorithm (Strebelle 2002). The workflow is generic and any other geostatistical sampler can be applied. Alternatively, we could simply ignore the interfaces mapped along the breakout intervals. Each type of well log is also attributed with a certain weight, which balances its reliability and weighs the amount of information it contains. The sum of all weights is equal to 1 and the individual weights are model parameters to be calibrated.

### Block 2: mapping interfaces

The mapping of interfaces is carried out by the algorithm of Lanning & Johnson (1983). Enhancements have been made to (1) speed up the calculations to enable on-site decisions and (2) to accommodate the calibration of the model parameters, reducing the subjectivity in the mapping process. The algorithm uses a set of Walsh functions to further discretize the signals transformed in block 1 (usually fuzzy and suffering from many small-scale fluctuations; see Fig. 2a, c). Walsh functions (Walsh 1923) are typically used in image rendering (Kennedy 1971 and references cited therein). They form an ordered set of rectangular ortho-normal waveforms and are described in detail in the original paper by Lanning & Johnson (1983). Low-pass signal filtering (here, a synonym to signal reconstruction) results in stepwise segmented versions of the transformed signals with enhanced, easy to visualize, transitions (solid red lines in Fig. 2c, d). The segment width is controlled by the number of Walsh functions and the chosen discretization. A low number of Walsh functions leads to an over-smoothed reconstruction of the signal (e.g. the first function represents the mean value; Fig. 2d) where trend changes cannot be identified. By contrast, a large number of Walsh functions may lead to a spiky reconstructed signal similar to the transformed signal, where identifying trends is difficult. The segment width is directly related to the mapping accuracy through the so-called minimum resolvable layer thickness (Lanning & Johnson 1983) and rock units thinner than this magnitude may not be resolved. The step width also controls the CPU time required to map interfaces. For all these reasons, the discretization and number of Walsh functions are model parameters to be calibrated.

Given a set of segmented signals with corresponding weights, the mapping algorithm proceeds as follows. Taking the first depth (in general, the length along the borehole) as a fictitious interface, the algorithm moves to the next segment along each log and computes the difference between the current standardized value and the mean on all previous segments since the last interface common to all logs. These individual differences (one per log) are then multiplied by the corresponding weight. If a certain log was not recorded at that depth, or if the measurement was corrupted, the corresponding individual difference is simply ignored. The sum of weighted differences yields the so-called pick value. Large pick values denote large changes occurring at all logs (on average if the weights are uniform). The pick value is then compared with another model parameter, the so-called check value. If the pick value is greater than, or equal to, the check value, then the current depth is picked as an interface and the algorithm restarts from the current depth. The check value is the most important model parameter because it controls the amount the picked facies (i.e. high check values lead to fewer picked interfaces). The check value is also the most uncertain parameter because it depends on the number and type of intervening logs. Logs with few trend changes require lower check values to identify a given number of interfaces. As such, the check value must be calibrated. Figure 3 shows the sensitivity of the mapping algorithm to the check value by comparing the interfaces picked using check values of *w* = 0.12 and *w* = 0.16 with those identified by a standard interpretation of the composite of logs and/or interface depths from lithostratigraphic profiles, termed here measurements. Two conclusions become apparent. First, the algorithm identifies many interfaces when the check value is low (a check value *w* = 0 would lead to the picking of an interface at each segment of the reconstructed logs). Second, the algorithm identifies most measured interfaces with accuracy (some of them within one-half of the segment width). In fact, most interfaces picked by the mapping algorithm, but not recognized by the standard interpretation, are clearly visible by closely examining the existing cores at the centimetre scale.

### Calibration of model parameters

The choice of model parameters plays an important part in the accuracy of the mapping. The discretization and number of Walsh functions control the degree of filtering (and correspondingly the amount of information retained from the original logs). The weights control the amount of information gathered from the reconstructed logs. The check value controls the amount and location of the picked interfaces. Provided that some measurements or benchmarks of the interface depths are available, the problem of finding the optimum values of model parameters can be defined as an inverse problem. The inverse problem or, in broad terms, inverse modelling, refers to the process of gathering information about the model from measurements of what is being modelled (Carrera *et al.* 2005). Measurements come from two different non-exclusive sources: (1) on-site, from the driller or the field geologist, and/or (2) from a posterior standard interpretation of the well logs (or a combination of both).

Model parameters can be calibrated manually or automatically. Automatic calibration frees the modeller of the burden of having to deal with a tedious manual trial and error process of fine-tuning the model parameters to achieve a best fit of the available measurements. Automatic calibration can be carried out by brute force (i.e. simply exploring different combinations of parameters) or by gradient search methods. In this work, the small number of model parameters allowed us to explore the parameter space. A mean square error function penalizing large differences between the model outputs and measurements is proposed:
*N _{m}* is the number of measurements

*m*,

_{i}*c*(

_{i}**) is the automatically picked interface closest to**

*p**m*, which depends on the model parameters (e.g. the check value or log weights) arranged in vector

_{i}**and**

*p**λ*are the weights characterizing the reliability of individual measurements (all set to 1 in this study). The vector

_{i}**contains the discretization**

*p**d*, the number of Walsh functions

*N*, the check value

*w*and the vector of weights balancing the information contained in each log. The second term in equation (1) penalizes the picking of a large number of interfaces

*N*, which would facilitate the calibration, and is weighed by a ponderation factor

_{c}*β*(set to 0 in this analysis). As such, very small values of the check value, leading to picked interfaces at almost each segment in the reconstructed signals, are also penalized.

A prior analysis revealed the low sensitivity of the mapping algorithm to both the discretization *d* and the number of Walsh functions *N*. The discretization must be dense enough to take full profit of the available log data. In the same line of arguments, the number of Walsh functions defining the segmented logs must be sufficiently high. In this analysis, and for ease of visualization, we set the discretization *d* to 25 cm (the measurement intervals were 1 and 10 cm for natural and spectral gamma ray logs, respectively) and the number of Walsh functions *N* to 304 (15% of the maximum, 2048). Tests made with *d* = 10 cm and *N* = 2048 (no segmentation of signals) did not lead to significantly different results. Figure 4 shows the sensitivity of the mapping algorithm to the check value and the weight of the natural gamma ray log *α _{N}* (the weight of the spectral gamma ray log is obtained directly by difference, 1 −

*α*). The contour curves of the objective function in equation (1) reveal a minimum using a check value of 0.1 and weights of 0.24 and 0.76 for the natural and spectral gamma ray logs, respectively. Figure 5 shows the interfaces mapped using these parameters. A visual comparison between Figures 3 and 5 shows the importance of calibration. In Figure 5, all the interfaces from the stratigraphic profile are detected by the algorithm with a mean error of 1.3 m (2.73 and 3.22 m for check values

_{N}*w*= 0.12 and

*w*= 0.16, respectively, in Fig. 3).

### Block 3: including facies type and hydrogeological properties in the calibration

The objective function in equation (1) can easily be enriched with additional terms of the same kind to account for other types of measurements or prior information about the model parameters. In this study, only depths of interfaces have been used. However, the workflow is generic and can accommodate other measurements, e.g. facies types (accounting for mineralogical contents), shale volume, porosity measurements or hydraulic conductivity values arising from laboratory experiments, separate interpretations of hydraulic tests or prior knowledge. Facies type can be detected automatically if the thorium and potassium contents defining the spectral gamma ray counts are provided (Fig. 6a). The calculation of shale volume *V*_{shl} involves a transformation of the spectral gamma ray log:
*SGR*, *SGR*_{cln} and *SGR*_{shl} are the spectral gamma ray values measured in the borehole, in a clean sand and in a shale zone, respectively. *SGR*_{cln} and *SGR*_{shl} are also model parameters to be calibrated, although a first estimate can be obtained from the cross-plot of the solid density and spectral gamma ray. Figure 6b shows the distribution of the clay content in the BDB-1 borehole attained with *SGR*_{cln} = 10 API and *SGR*_{shl} = 90 API (not obtained from actual well logs). In addition, the porosity can be calculated from certain well logs (e.g. density or combinable magnetic resonance) via empirical relationships and additional model parameters. In a second step, the hydraulic conductivity can be inferred from porosity through, e.g. the Kozeny–Carman law (Kozeny 1927). These distributions can then be related to prior knowledge and/or actual laboratory or field measurements to further constrain the model parameters. Secondary information on the facies type, mineralogy or the relevant hydrogeological properties enriches the accuracy of the mapping and enhances the characterization of the stratigraphic sequence. Unfortunately, considering such measurements involves the calibration of additional parameters, making the brute search calibration process slower (although still feasible). In such a case, gradient search methods (e.g. the regularized pilot points method; Alcolea *et al.* 2006) should be used. However, this is beyond the scope of this paper. Instead, we focus on the parameters defining the accuracy of interface mapping.

### Block 4: geostatistical post-processing

Once the distributions of the shale volume, porosity and hydraulic conductivity have been estimated, variograms describing the variability of such properties can be inferred, which allow us to establish the correlations between different boreholes and enables geostatistical modelling (Fig. 6c, d). Variograms are obtained automatically by means of jack-knife techniques (Quenouille 1949).

The workflow can also be used to address uncertainties in the estimated parameters. These arise from the propagation of the inherent uncertainties of well logs (e.g. caused by environmental effects, mud density or other statistical effects related to natural decay) or of the empirical functions used in the transformations. To that end, the workflow can be cast in the framework of Monte Carlo statistical analysis (Yang & Torres-Verdín 2015).

## Applications

This section presents two applications involving different scales and dimensionalities. The first application consists of the one-dimensional mapping of lithostratigraphic units at the metre to decametre scale in 13 boreholes in northern Switzerland. The second application consists of the characterization of the heterogeneity of a small rock surface at a wall in the so-called TT niche of the Mont Terri URL. Thus the analysis is two-dimensional (although made as a composite of one-dimensional interpretations) and at the millimetre to centimetre scale. The objective is two-fold. First, we show that the workflow is capable of mapping facies with accuracy regardless of the scale. Second, we test the sensitivity of the mapping algorithm to its main parameter, i.e. the check value.

### Mapping of large-scale stratigraphic units

The location of the 13 boreholes under study, with depths varying in the range 150–2500 m, is shown in Figure 7. Several logs were collected at each borehole and the composite was interpreted to compare the geostatistical results with the interfaces obtained from lithostratigraphic profiles (lithostratigraphic profile from Lindau 1 published in Büchi *et al.* 1965; from Lausen in Vogt *et al.* 2016; from Benken in Nagra 2001; all other values were taken from stratigraphic profiles in Nagra 2014; all geophysical data was obtained from Nagra archives). Interfaces from stratigraphic profiles are used here as measurements to calibrate the check value. For simplicity, one single representative log has been chosen at each borehole. The mapping of interfaces is carried out twice. First, a common check value of *w* = 0.16 is set (i.e. calibration is not performed). This high value has been purposefully chosen to illustrate the impact of a wrong choice of the check value.

Table 1 shows the calibrated values and mean mapping errors, defined as the absolute differences between the measured interfaces and the closest calculated interface depth. The table also reports the number of interpreted interfaces not detected by the mapping algorithm. To that end, a small threshold error of 1 m was considered. The first conclusion that becomes apparent is the negative impact of a wrong choice of the check value on both the mapping accuracy and the number of missed interfaces. If a wrong (in the sense of not calibrated) check value of *w* = 0.16 is used, then the algorithm misses 36% of the total interfaces to be mapped (344 overall in all 13 boreholes). Correspondingly, the mean mapping errors are high. These values decrease dramatically by calibrating the check value (2% of missed interfaces and small individual mapping errors). In this instance, the mean error (all boreholes) is 0.6 m only. This is especially relevant bearing in mind how the measurements are obtained (i.e. by visual and subjective inspection of a composite of logs) and because only one representative log was used in the calibration. The strong variability of the calibrated check values in the range 0.02–0.15 is explained by local heterogeneities and by the different accuracies of the logs used in the calibration.

Figure 8 shows a comparison between the gamma ray measurements of acoustic borehole image (ABI) logs in boreholes Hemmental-2 and Tegerfelden-2, along with the segmented signals and mapped/measured interfaces. The ABI log measured in Hemmental-2 shows many more small frequency fluctuations than that in Tegerfelden-2. These fluctuations are alleviated, but do not completely disappear, after segmenting the signal (304 Walsh functions were used in all instances). Because the algorithm aims to pick an interface close to the measured interfaces, it requires a lower check value when the amount of fluctuations in the segmented signal is small. Figure 9 shows a box-plot of mapping errors for all boreholes. The mean mapping error is <1 m in all instances, except in Schafisheim (1.15 m). These errors are very small compared with the mean facies thickness (Table 1).

### Mapping of small-scale features

A hand-held spectrometer monitored the spectral gamma emission along 14 horizontal scanlines over a *c.* 75 × 15 cm^{2} surface at a rock wall in the TT niche of the Mont Terri URL (Fig. 10). Several structures (darker colours, dipping 45°) oriented parallel to the bedding planes can easily be identified. A closer look at these structures (Fig. 11) shows both their elongated shape (thickness <10 mm in most instances) and a marked intra-facies variability, defined by interbedded microstructures mixing different mineralogies. The manual lifting of the spectrometer using a forklift made it difficult to properly georeference the horizontal scanlines with respect to the photographs in Figures 10 and 11, making it impossible to calibrate the model parameters with the image. Instead, a blind prediction of the shape and distribution of the structures is inferred from spectral gamma ray readings.

First, the mapping algorithm is separately applied to each scanline. After a first sensitivity analysis, we use a check value *w* = 0.02, a discretization *d* = 1 cm corresponding to the horizontal and vertical measurement frequency, and *N* = 127 Walsh functions to segment the spectral gamma logs. Overall, 247 interfaces are mapped. For each identified facies, the mean thorium to potassium ratio is calculated and categorized (Fig. 12). This variable is considered as a good indicator of the facies type. Visual inspection of Figure 12 shows certain structures dipping at 45°. However, the small thickness of the measured features cannot be captured by the mapping algorithm because the discretization is coarse and the spectrometer averages the gamma readings across the scintillator crystal (i.e. a surface of 7.5 × 7.5 cm^{2}).

In an attempt to reduce the impact of the measurement frequency (1 × 1 cm^{2}), the five thorium to potassium ratio categories in Figure 12 are used as conditioning measurements to estimate the distribution of facies. Variograms are first estimated in the directions of 0, 45 and 90° (Fig. 13). In view of the preferential alignment of structures parallel to bedding planes dipping at 45°, it is not strange that the best variogram fit is attained along that direction. In fact, experimental variograms along horizontal and vertical directions are fuzzy and characteristic of the variability of a non-stationary variable (Samper & Carrera 1996). The best-fitting model variogram is spherical, with a nugget 0.05 and a sill 0.85. The non-zero nugget confirms the presence of scales of variability smaller than the distance between measurements (10 mm). Such scales are easily observable in Figure 11. The longitudinal range and the variogram bandwidth are 90 and 10 mm, respectively, coherent with the mean length and thickness of some structures. The domain is now discretized using a regular grid with a step size of 1 mm and kriging of the mapped categories in Figure 12 is carried out using the model variogram in Figure 13. The results, in terms of the estimated categories and corresponding kriging variances, are shown in Figure 14. As observed, the structures are correctly aligned with the bedding planes. Obviously, this is a byproduct of the use of the 45° oriented variogram. The identified structures resemble the shape of those in Figures 10 and 11. Unfortunately, the identified structures cannot be unambiguously correlated with those in the photographs due to the lack of a common georeference system. However, global trends in the alignments are clearly visible in the estimated field. This shows that the workflow is capable of identifying small-scale features if these are properly captured by the measured logs.

## Conclusions

A workflow for the automatic interpretation of well logs (or, in more general terms, measurements with one correlating variable) has been presented. It integrates a mapping algorithm, inspired by that of Lanning & Johnson (1983), and empirical functions to identify lithostratigraphic units and derive hydrogeological properties, enhancing the standard interpretation of stratigraphic sequences in boreholes. Efforts have been devoted to: (1) the calibration of model parameters, making the interpretation of well logs objective and traceable; (2) to the algorithmic efficiency, facilitating on-site decisions based on actual data; and (3) to the inclusion of geostatistical tools for post-processing, enabling further modelling exercises, e.g. correlations between boreholes.

The calibration exercises presented here are illustrative and simple and focus on the check value, which controls the accuracy of the mapping, and on the weights balancing the amount of information contained in each well log. However, the workflow is generic and allows us to estimate the remaining parameters (optimum discretization and number of Walsh functions used in signal enhancement, and possibly others governing empirical functions for log transformation). In the same line of argument, we have only used the depth of facies interfaces obtained by a standard interpretation of the composite of logs as measurements. Interfaces detected during drilling and/or interpretations of laboratory or hydraulic tests (e.g. shale volume, facies type, porosity and/or hydraulic conductivity) can also be used as measurements to further constrain the model parameters and enrich the description of the stratigraphic sequence.

The suggested workflow was applied to two different problems involving different scales and dimensionalities. The results of the first application show the importance of calibration. Optimum model parameters lead to small mapping errors and therefore to more accurate descriptions of the large-scale stratigraphic sequence in 13 boreholes in northern Switzerland. After calibration, the workflow identifies with accuracy 98% of the large-scale facies interfaces present in the 13 boreholes (only 7 missed out of 344 measured interfaces; Table 1). Instead, it identifies only 64% with wrong parameters (123 missed). The versatility of the workflow to resolve different scales is addressed by the second application. The workflow was used to map facies along horizontal scanlines over a small surface containing small-scale structures with thicknesses <1 cm. The absence of a georeference system precluded the calibration of the model parameters, impeding the resolution of thin structures. The discretization issue was addressed by estimating a kriged field from the facies type identified by the workflow. Although the reconstructed field could not be directly compared with on-site observations, the alignments and shapes of the identified structures resemble those observed in the field.

Much remains to be done. The model parameters are calibrated by exploring the parameter space (i.e. the brute force testing of combinations of model parameters), which has a negative impact on CPU time. This issue can be tackled by using gradient search methods. The impact of other type of logs or measurements (e.g. facies type or shale volume) must be further explored to validate the use of empirical functions in the workflow. Uncertainties in well logs (e.g. those caused by environmental effects, mud density or statistical uncertainties due to natural decay) can be taken into account by casting the workflow in a Monte Carlo simulation framework. The results presented here should be viewed as a hopeful step in the direction of automatically and efficiently interpreting geophysical well logs.

## Acknowledgements

The first author developed part of this work during his time with TK Consult AG. The valuable review comments and suggestions by Thomas Spillmann (Nagra), Jorge Jodar, A. Soler (University of Barcelona, Spain) and J. Jimenez-Martinez (ETHZ-EAWAG, Switzerland) to improve the technical quality of this paper are highly appreciated.

## Funding

This work has been supported by the Swiss National Cooperative for the Disposal of Radioactive Waste (Nagra) and the Mont Terri Consortium.

- © 2018 The Author(s). Published by The Geological Society of London. All rights reserved