Using a multiple variogram approach to improve the accuracy of subsurface geological models 1

Subsurface geological models are often used to visualize and analyze the nature, geometry, and variability of geologic and hydrogeologic units in the context of groundwater resource studies. The development of three-dimensional (3D) subsurface geological models covering increasingly larger model domains has steadily increased in recent years, in step with the rapid development of computing technology and software, and the increasing need to understand and manage groundwater resources at the regional scale. The models are then used by decision makers to guide activities and policies related to source water protection, well field development, and industrial or agricultural water use. It is important to ensure that the modelling techniques and procedures are able to accurately delineate and characterize the heterogeneity of the various geological environments included within the regional model domain. The purpose of this study is to examine if 3D stratigraphic models covering complex Quaternary deposits can be improved by splitting the regional model into multiple submodels based on the degree of variability observed between surrounding data points and informed by expert geological knowledge of the geological– depositional framework. This is demonstrated using subsurface data from the Paris Moraine area near Guelph in southern Ontario. The variogram models produced for each submodel region were able to better characterize the data variability, resulting in a more geologically realistic interpolation of the entire model domain as demonstrated by the comparison of the model output with preexisting maps of surficial geology and bedrock topography as well as depositional models for these complex glacial environments. Importantly, comparison between model outputs reveals significant differences in the resulting subsurface stratigraphy, complexity, and variability, which would in turn impact groundwater flow model predictions. Résumé : Des modèles de géologie de subsurface sont souvent utilisés pour visualiser et analyser la nature, la géométrie et la variabilité d’unités géologiques et hydrogéologiques dans le contexte d’études sur les ressources d’eau souterraine. Le développement de modèles géologiques de subsurface en trois dimensions (3D) couvrant des domaines modélisés de plus en plus grands s’est accéléré ces dernières années, parallèlement à l’essor rapide de la technologie informatique et des logiciels et au besoin accru de comprendre et de gérer les ressources d’eau souterraine à l’échelle régionale. Les modèles sont ensuite utilisés par les décideurs pour orienter les activités et politiques touchant à la protection des sources d’eau, la mise en valeur de champs de captage et l’utilisation industrielle ou agricole de l’eau. Il importe de s’assurer que les techniques et procédures de modélisation puissent délimiter et caractériser de manière exacte l’hétérogénéité des divers milieux géologiques contenus dans le domaine régional modélisé. Le but de l’étude est d’évaluer si les modèles stratigraphiques 3D couvrant des dépôts quaternaires complexes peuvent être améliorés en divisant le modèle régional en plusieurs sous-modèles en fonction du degré de variabilité observé entre les données d’un secteur et à la lumière de connaissances géologiques de pointe sur le cadre géologique– sédimentaire. Une démonstration est faite en utilisant des données de subsurface de la région de la moraine de Paris, près de Guelph, dans le sud de l’Ontario. Les modèles de variogramme produits pour chaque région de sous-modèle sont en mesure de mieux caractériser la variabilité des données, ce qui se traduit par une interpolation plus réaliste sur le plan géologique de l’entièreté du domaine modélisé, comme le démontre une comparaison des sorties du modèle à des cartes préexistantes de la géologie de subsurface et du relief du substrat rocheux, ainsi que des modèles sédimentaires pour ces milieux glaciaires complexes. Fait à noter, la comparaison de sorties de modèle révèle des différences significatives dans la stratigraphie, la complexité et la variabilité qui en découlent, ce qui devrait avoir une incidence sur les prédictions de modèles d’écoulement de l’eau souterraine. [Traduit par la Rédaction] Received 26 June 2016. Accepted 29 July 2017. K. MacCormack* and E. Arnaud.† School of Environmental Sciences, University of Guelph, Guelph, ON N1G 2W1, Canada; G360 Institute for Groundwater Research, University of Guelph, 360 College Avenue, Guelph, ON N1G 2W1, Canada. B.L. Parker. G360 Institute for Groundwater Research, University of Guelph, 360 College Avenue, Guelph, ON N1G 2W1, Canada; School of Engineering, University of Guelph, Guelph, ON N1G 2W1, Canada. Corresponding author: Kelsey MacCormack (email: kelsey.maccormack@aer.ca). *Present address: Alberta Geological Survey, Alberta Energy Regulator, 4999-98 Avenue, Edmonton, AB T6B 2X3. †Emmanuelle Arnaud currently serves as a Special Editor; peer review and editorial decisions regarding this manuscript were handled by Special Editor


Introduction
Many three-dimensional (3D) subsurface geological models were created for various regions in southern Ontario following the Walkerton Inquiry (O'Connor 2002) and the consequent source water protection initiatives.Such models have been used in Ontario and elsewhere for a number of environmental management applications, including risk assessments of municipal or hazardous waste and the management and protection of groundwater resources (e.g., Artimo et al. 2004;Bajc et al. 2011;Thomason and Keefer 2011;Atkinson and Glombick 2015).In other cases, 3D geological models have been used in the context of the delineation and remediation of contaminants, and the exploration and extraction of aggregate and mineral resources (Houlding 2000).The overall aim of these models is to characterize the nature and distribution of subsurface materials, which typically necessitate modelling near-surface glacial deposits, to constrain the modelling of the groundwater flow system, contaminant migration pathways, transport rate of contaminants, and risk to receptors, or to optimize natural resource extraction (Wycisk et al. 2009;Babak et al. 2014;Dunkle et al. 2016).Accurately characterizing glacial deposits can be particularly difficult, as subsurface stratigraphy in glacial settings typically have short length scales, abrupt facies transitions, and erosionally truncated relationships as a result of the dynamic nature of these depositional environments and multiple and spatially variable erosional and depositional events associated with repeated ice advances and retreats.
The scale of 3D subsurface geological models developed in the context of groundwater resource management has been rapidly increasing over the past decade to cover larger model domains (Bajc and Newton 2005;Kessler et al. 2005;Logan et al. 2006;Arihood 2008;Gunnink et al. 2013;McLaughlin et al. 2015).With increasing size, regional models (hundreds of square kilometres or more) are likely to encompass a variety of depositional environments, resulting in multiple heterogeneous zones within a stratigraphic-hydrostratigraphic model (Weissmann et al. 1999;Babak et al. 2014;Dunkle et al. 2016).Modelling heterogeneity across regional study areas can lead to erroneous interpretations when models are constructed using geostatistical interpolators that attempt to incorporate this variability without acknowledgment of these subdomains (de Marsily et al. 2005).
The most commonly applied geostatistical interpolation algorithm is ordinary kriging (Kravchenko and Bullock 1999;Johnston et al. 2001;Jones et al. 2003;Kravchenko 2003;Mueller et al. 2004).Ordinary kriging is a distance weighted estimation algorithm that uses a variogram to characterize the variability between individual data points across a study area to optimize the weights assigned to the data points during the estimation of each grid node (Isaaks and Srivastava 1989;Krajewski and Gibbs 1996).If there is a large amount of spatial variability within the dataset, it may not be possible to fit a reliable variogram, which will result in poor estimates with large associated errors (Weber and Englund 1994;Gringarten and Deutsch 2001).The problem with interpolating subsurface units across regional models containing zones of both high and low variability (the current standard procedure) is that the resulting variogram must characterize the variability within the entire model area, and as a result, produces a blended estimate with unknown bias of the actual variability (Houlding 2000).This problem is compounded by the fact that geostatistical interpolators are often utilized as a "black box" tool by users who do not fully understand how the predictions were calculated or if they are realistic in their geologic context.Together, these factors may severely compromise the accuracy of the model outputs (Goodchild and Haining 2004).Unfortunately, little attention (if any at all) is paid to variogram analysis, which lies at the center of all geostatistical analysis and estimation for assessing how similar the existing and interpolated data are to one another (Gringarten and Deutsch 2001;Chiles and Delfiner 2012).
The objective of this study is to determine if a model can be improved by splitting the regional model into multiple submodels based on the degree of variability observed between surrounding data points and our understanding of the depositional environment.This modelling approach will aim to incorporate the existing geologic framework ("hard" data, e.g., topography, geomorphic features, geologic logs) with "soft" qualitative data (e.g., geological knowledge of specific depositional environments) to define the submodel domains, thus optimizing the measured data values as input to the geostatistical model.The hypothesis is that by separating the regional model domain into submodels based on the data variability and our understanding of the geologic context, a variogram can be produced for each submodel region that is more representative of the true variability within the data for each area, thus resulting in a more geologically realistic interpolation and model results for both the subdomains and entire model domain (Weissmann et al. 1999;de Marsily et al. 2005).
The regional model domain for this study encompasses a 165 km 2 area that includes a section of the Paris Moraine near Guelph, Ontario (Fig. 1).This was an ideal dataset to use for this study because the model area contains several geological environments and has been identified by the province of Ontario as an area requiring a water budget evaluation and development of a source water protection policy (Ontario Ministry of the Environment 2009).Developing methods to more accurately delineate the form and geometry of the subsurface sedimentary deposits as presented here can be used to help quantify the impact of the Paris Moraine on groundwater recharge and flow variability as well as inform our understanding of potential contaminant pathways and receptorsboth important aspects of current water resource management in the region (Ontario Ministry of the Environment 2009; Russell et al. 2009).The method developed here can also be used elsewhere, in other glaciated regions where a variety of depositional environments are expected to create significant heterogeneity across a modelling domain.

Background geology
The Guelph area was glaciated multiple times during the Quaternary period, resulting in a complex interfingering of glacial deposits overlying Paleozoic bedrock and multiple glacial landforms at surface (Karrow 1963(Karrow , 1968(Karrow , 1987;;Chapman and Putnam 1984;Barnett 1992).Drift thickness over bedrock ranges from 10 to 15 m in the outwash and till plains and up to 40 m along the Paris Moraine.The oldest glacial deposits in the region are defined as Early to mid-Wisconsinan in age and are either till (Canning Till) or stratified sediments, though these are only documented in a few places.Overlying these older deposits is a relatively continuous regional basal till (Catfish Creek Till) deposited during the major ice advance of the Nissouri Phase (Late Wisconsinan; Table 1) has been described by Karrow (1968Karrow ( , 1987)).The Catfish Creek Till unit is overlain by kame and outwash deposits of the subsequent Erie glacial retreat phase.The Port Stanley Till is found at surface in the NW section of the study area and is thought to be associated with the younger Port Bruce ice advance, whereas the Wentworth Till at surface in the SE sector records an ice readvance during the overall retreat of the Erie-Ontario lobe (Mackinaw Phase) (Table 1; Karrow 1987;Barnett 1992;Karrow et al. 2000).In some areas, Karrow (1968Karrow ( , 1987) ) described an additional till between the Catfish Creek Till and the tills at surface, which he interpreted as Port Bruce in age (Maryhill-Middle till), although the similarities between Maryhill and Port Stanley tills in this area, and the limited stratigraphic control and radiometric data, do not always allow differentiation of tills at depth.Lastly, laterally discontinuous kame and outwash deposits occur at surface throughout the study area and record glaciofluvial processes associated with ice retreat from the region.The most prominent landform at surface is the Paris Moraine, a NE-SW-trending topographic high that cuts across the study area (Sadura et al. 2006) and records stagnation or minor readvance during the Mackinaw Phase (Karrow 1968).It is defined by a core belt of hummocky topography and is composed of silty sandy diamict (Wentworth Till), as well as coarse-grained outwash, kame, and ablation deposits.The surrounding landscape consists primarily of outwash plains on either side of the moraine, as well as a drumlinized till plain with Port Stanley Till at surface on its NW side (Karrow 1968;Chapman and Putnam 1984).The Paris Moraine is thought to be an important hydrogeologic feature for enhanced recharge to the underlying bedrock aquifer and provides spatially distributed baseflow to adjacent creeks (Blackport et al. 2009;Russell et al. 2009).The Paris Moraine is also considered to be an important aggregate resource (Ontario Ministry of the Environment 2009).Delineating the 3D spatial geometry of the Paris Moraine and associated deposits is important for accurately assessing its economic and hydrogeologic value and potential, which in turn can be used to inform source water protection and groundwater resource management decisions.

Data compilation
A total of 1508 data points recording depth and occurrence of seven different stratigraphic units were integrated to produce this 165 km 2 model.Two hundred and sixty-three of these data points came from moderate-to high-resolution geologic logs from academic studies (8), and geotechnical reports (255), whereas the remaining 1245 boreholes were obtained from the Ontario Ministry of the Environment water well database (Fig. 1).This resulted in a data density of 9.1 data points / km 2 for the entire study area.When using lithology data to interpolate stratigraphic models, it is best to use only (or primarily) high-quality and reliable data.However, in this and many other cases, the use of low-quality data from more regionally extensive databases is necessary to improve the distribution of data across the study area, as there were very few high-quality data points, and many of the moderate-quality data were typically very shallow (7 m or less) and clustered around local construction sites.As a result, the more regionally extensive water well database (low-quality data) was used to provide information about the subsurface lithology where high-quality data were not available (MacCormack and Eyles 2010).To help reduce the number of inaccurate water well data included within the model, only data with a location accuracy of less than 100 m, and elevation accuracies of approximately 3 m or less were included.This methodology of filtering water well data has been applied in other studies forced to rely heavily on water well data (Gao et al. 2006;Puckering 2011).Removing inaccurate water well data prior to model interpolation has been done in previous studies such as Logan et al. (2006).

Lithologic and stratigraphic coding
The sediment description data provided for each borehole record were reduced from 108 unique descriptions to 22 textural categories.Combining similar sediment descriptions (such as sand, fine sand, sandy, silty sand, medium sand) into categories is necessary to gain a better sense of the significant textural variability across a study area (MacCormack et al. 2005;Dumedah and Schuurman 2008;Dunkle et al. 2016).In this study, 22 textural categories were then grouped into assemblages representing seven stratigraphic units: bedrock, undifferentiated tills, Port Stanley Till, ablation deposits, Wentworth Till, kame deposits, and outwash deposits.This was done by comparing the descriptions within the boreholes to the regional stratigraphic framework and assigning stratigraphic units depending on the thickness and sequence of sediment descriptions within each borehole (Fig. 2; Table 1).Bedrock was assigned to the base of boreholes that indicated bedrock or solid rock had been encountered and was drilled for at least 1 m.The undifferentiated till unit was used where till or diamict identified in the borehole records did not fit stratigraphically or lithologically with the other till units.Karrow (1968Karrow ( , 1974) ) suggests that finding till deposits older than the Port Stanley or Wentworth tills within the Paris Moraine is likely.However, the water well data typically do not contain sufficiently detailed sediment description to accurately identify and differentiate these older till or diamict units within the existing regional till stratigraphic framework (Table 1).The Port Stanley Till is identified where the borehole records identified a silty sandy till, or similar description at surface or at a relatively shallow depth below surface and in the NW sector of the study area (Table 1; Karrow 1968;Russell et al. 2009).Ablation deposits were identified as stratified sand and gravel typically found in a stratigraphically similar position as the Port Stanley Till (Table 1).The Wentworth Till is a sandy to silty till unit and is typically found in the SE sector of the study area (Table 1; Karrow 1968;Russell et al. 2009).Above the Wentworth Till and primarily along the moraine and back slope in the SE sector are numerous isolated and discontinuous packages of sand and gravel that have been identified as kame deposits (Table 1; Cowan 1975;Karrow 1987).The uppermost stratigraphic unit is a more extensive package of stratified sand and gravel identified throughout the study area as outwash deposits associated with the stagnation and final retreat of the Laurentide Ice Sheet from this region (Table 1; Karrow 1987;Russell et al. 2009).All borehole data were coded with picks identifying the elevation of the top contact for each stratigraphic unit (Fig. 2).If a stratigraphic unit was determined to be missing within a borehole, the unit's top and base elevations were assigned the top elevation of the underlying horizon, resulting in a zero thickness, except if the missing unit(s) occurred at the bottom of the well, in which case they were not recorded.

Kriging and the use of variography
Borehole data were interpolated using the kriging algorithm in RockWorks15 (a geological modelling software package; RockWare Inc. 2009) to predict the vertical and horizontal distribution and overall geometry of each stratigraphic unit.Owing to the substantial differences in the spatial extent of the stratigraphic units, it was necessary to produce modelled surfaces for the top and base of each unit.For example, the kame deposit is often laterally discontinuous so that its top surface cannot always be used to define the base of the overlying outwash deposit.
Kriging is the most common algorithm used to interpolate geologic data (Kravchenko and Bullock 1999;Johnston et al. 2001;Jones et al. 2003;Kravchenko 2003;Mueller et al. 2004).This geostatistical algorithm is capable of considering the spatial proximity of the observed data values (i.e., depth and occurrence of the seven stratigraphic units) as well as incorporate statistical properties of the measured observations into the calculations of the predicted value (Isaaks and Srivastava 1989).This is achieved by Table 1.Generalized stratigraphy as described by previous studies in this region.utilizing a variogram to provide additional information on the spatial autocorrelation of the observed stratigraphic units with respect to the model grid nodes (Isaaks and Srivastava 1989;Johnston et al. 2001).This provides an effective method of quantifying the similarity (or dissimilarity) between data points with respect to distance (Cressie 1993;Logan et al. 2006;Chiles and Delfiner 2012).
The variogram is defined by the following equation: The basic function of a variogram is to illustrate the semivariance between two data points (X at locations {x i } and {x i+h }) over varying lag distances (h) (Houlding 2000).This equation calculates the sum of the squared differences between pairs of points separated by distance.The number of points used in the comparison is n; therefore, n − h represents the number of comparisons (Davis 2002).
Variogram components that provide information about the spatial autocorrelation and continuity of the data (such as nugget, major range, direction of major range, sill, type of variogram model, and correlation coefficient) were analyzed to assess the impact of subdividing the regional model on the interpolation of each stratigraphic unit within each submodel.

Subdividing the model
Interpolating and modelling complex geological environments within a single model is often very difficult, as such deposits typically include a substantial amount of natural (real) variability that modelling algorithms are unable to account for (Logan et al. 2006;Keefer 2007).Conventional geostatistical algorithms, such as kriging, require the assumption that the data being interpolated is stationary, such that the variance and mean are the same throughout the model domain (Weissmann and Fogg 1999;Davis 2002).When interpolating data across a glaciated region, it is unlikely that the assumption of stationarity will remain valid considering the dynamic nature of these depositional environments that result in both variable sediments at surface and a complex subsurface stratigraphy (Weissmann and Fogg 1999;Babak et al. 2014).Thus, to improve the results of the regional model interpolation, the study area was subdivided into regions of high and low complexity (variability) that were analyzed independently.The regional model area was divided into submodels (representing zones of local stationarity) by assessing the data variability and geological complexity within the context of the regional conceptual model.The original regional model domain was subdivided into three submodel zones (Fig. 1) by comparing changes in the lateral facies variability at depth in the borehole dataset (Fig. 2) with stratigraphic and physiographic information from previous scientific publications (Table 1; Karrow 1968;Sadura et al. 2006;Russell et al. 2009).These publications typically include details on the spatial extent of these different environments at surface as well as typical stratigraphic units associated with each environment based on surficial mapping, borehole records, and ground-penetrating radar datasets.
Subsurface sediments in zone 1 are composed of a series of diamict units overlain in many areas by channelized fans of outwash sand and gravel.Zone 1, which encompasses outwash and till plain physiographic elements at surface, was identified as a submodel because the stratigraphic layers in this region tend to be more continuous and laterally extensive in comparison to the neighbouring hummocky deposits of the moraine.Zone 2 is typically very hummocky at surface with a complex internal stratigraphy of Wentworth Till at surface, older till(s) at depth, as well as isolated outwash, kame, and ablation deposits.Zone 2 broadly corresponds to the elongated form and geometry of the Paris Moraine and will likely result in data from the stratigraphic units being correlated in a direction along the length of the moraine and not in the direction of ice movement as expected for the outwash and till plain sediments (Fig. 1).Zone 3 encompasses the area to the southeast of the Paris Moraine.These deposits are distinct from the main moraine body and modelled independently owing to the highly variable and complex nature of the subsurface sediments and stratigraphic units in this area.The com- Pagination not final (cite DOI) / Pagination provisoire (citer le DOI) plexity of the zone 3 deposits (Fig. 1) within this region are the result of sediment reworking by gravity and glaciofluvial processes and drainage water ponding between the topographic highs of the Paris Moraine and nearby Galt Moraine.
Once the submodel domains were identified, polygon filters were created to subdivide the data based on the defined submodel areas (Fig. 1).The individual submodels were interpolated using the same process as the original regional model.The ordinary kriging algorithm was used to interpolate each model, and the variogram results were analyzed to see how the nugget, distance of major range (herein referred to as range), direction of range, sill, type of variogram model, and correlation coefficient varied for each submodel.A grid spacing of 100 m 32 × 100 m in the X and Y directions, respectively, was used for each model.This grid node spacing was chosen to ensure that our models were able to capture the smaller-scale structures of the kame and outwash deposits, while not being too computationally intensive.This produced a total of 615 600 grid nodes covering the 165 km 2 study area, which offered good resolution for defining the spatial extent of each stratigraphic unit.

Recombining the submodels
Once all the submodels were interpolated, the upper and lower grid surfaces of each stratigraphic unit were joined together using the Grid Editor in RockWorks15 (RockWare Inc. 2009).The stratigraphic units were joined with the respective stratigraphic units in neighbouring submodel domains to form continuous layers over the entire original regional model domain.There were generally sufficient data points along the common boundaries between the submodels to permit adjacent models to be seamlessly recombined.Where stratigraphic unit grids did not align across submodel boundaries, a local grid filter (4 grid cells × 4 grid cells) was applied to smooth the transition between the grids.Once the grid surfaces from neighbouring submodel areas were joined together, the regional model was rebuilt using the new grids for each stratigraphic unit (Fig. 3).

Comparing the model output results
To assess the impact that subdividing the regional model had on the model results, six variogram parameters were evaluated for the regional model and each of the submodels: (1) distance of range, (2) nugget, (3) sill, (4) type of variogram model, (5) direction of range, and (6) correlation coefficient of the variogram model.By subdividing the regional model into smaller zones within which the data are more similar, it is expected that the variogram produced will be more characteristic of the data and better able to model the geological features within each submodel.This is in turn should be reflected in changes in variogram parameters as discussed later in the text.The scale of the grids used to model the regional model and submodel are the same, and the density of data for each zone are quite similar (10.8, 7.8, and 8.5 data points / km 2 for zones 1, 2, and 3, respectively).Therefore, the change in the size of the modelled area in and of itself is not expected to affect the variogram parameters.
Root mean square error (RMSE) was used to compare the model output results from the standard original regional model and the proposed submodels with the original observed or measured data points.The RMSE is a good comparative statistic for assessing model output, as it provides a global indication of how similar the interpolated values are to the observed or measured data point values (MacCormack et al. 2013).When analyzing the RMSE statistics, a small RMSE value indicates that the interpolated values for the output model are more similar to the observed data point values, whereas a large RMSE value suggests that the interpolated model values are less similar to the observed data points.Thus, RMSE values are used here to determine how well the model fits the observed data values, with low RMSE values indicting a high degree of model accuracy (Zimmerman et al. 1999;Davis 2002;Dille et al. 2003;Jones et al. 2003;Mueller et al. 2004).
where Z ˆ(s i ) is the interpolated value at the point (s i ), z(s i ) is the observed (true) value from the original dataset at that same location, and n is the number of points within the input dataset.

Results and discussion
The goal of this study was to investigate if by subdividing the regional model domain, based on the data variability and our understanding of the geological framework and depositional environment, the variogram produced for each submodel would be more representative of the variability within the data for each area, thus resulting in a more geologically realistic interpolation of the entire model domain.The regional 3D stratigraphic reconstructions were modelled using the same algorithm for both the standard approach of modelling the entire regional domain at one time (original regional model) and the proposed methodology of combining submodels of the regional domain (subdivided regional model); these reconstructions were then compared with one another (Fig. 3).The kriging algorithm used to interpolate data in each of the models is highly sensitive to the variogram model parameters (Myers 1997).Therefore, the variogram parameters (distance of major range, nugget, sill, type of variogram model, direction of major range, and correlation coefficient) from the original regional model and each of the submodels are compared to provide additional insight into the difference in model output.

Comparing the nature and geometry of stratigraphic units
Visually, there are substantial differences between the two models.The subdivided regional model was able to characterize additional complexity and variability both at surface (Fig. 3) and at depth (Figs. 4,5).The surficial outwash deposits of the subdivided regional model are much less extensive, and constrained to the NW and SE corner of the model domain, which is more consistent with published surficial geology maps of the area (Fig. 3C; Ontario Geological Survey 2003).The fence diagram shows that other stratigraphic units at depth are also less extensive across the domain, resulting in more complex geometry of units that pinch out or thicken in specific areas and at specific depths within the domain (e.g., NW and EW panels in Fig. 4); this more complex geometry is thought to be more geologically realistic than the simpler "layer-cake" model of the original regional model, considering the sediment heterogeneity encountered in glacial deposits in general (Stephenson et al. 1988) and in this area in particular (Blackport et al. 2009;Russell et al. 2013;Arnaud et al. 2017).The variable geometry of these units is an important consideration in the context of both local and regional hydrogeological investigations where different unit geometry (with distinct lithology characteristics) can affect overall recharge rates of groundwater, travel time, and distribution of contaminants.Improving model predictions of the form and geometry of highly heterogeneous areas like this may help to identify where windows in aquitard materials exist, which may compromise the integrity of the aquitard and water quality in the underlying aquifer.
Notably, there are also significant differences in the geometry of the sediment-bedrock interface between the original regional model and the subdivided regional model outputs (Fig. 4, 5), with significantly more relief on the bedrock surface and better constrained buried bedrock valleys in the subdivided regional model compared with the original regional model.This variable sediment bedrock interface is consistent with the variable bedrock Fig. 3. Three-dimensional geological models of (A) the original regional model for which all the data were interpolated and (B) the subdivided regional model for which the three submodels were interpolated independently and then recombined.Models have a vertical exaggeration of 10×.(C) Surficial geology map of the study area for comparison with model outputs; note the improved correlation between the surficial geology map and the subdivided regional model.Note that orange and yellow in (A) and (B) represent the sand and gravel (undifferentiated) of the outwash and kame stratigraphic units, whereas orange and yellow in (C) represent gravel and sand, respectively.Map in (C) created using ESRI ArcMap 10.1 and data from Ontario Geological Survey (2003); white lines represent the submodel boundaries defined in Fig. 1; black hatched lines represent the location of the cross sections in Fig. 4. [Colour online.]surface topography and the relatively deep and well-defined buried bedrock valleys documented previously (Karrow 1968(Karrow , 1987;;Gao et al. 2006;Cole et al. 2009) and observed in recently drilled boreholes within the study area (Steelman et al. 2017).Relief on the bedrock topography surface is important in both local and regional groundwater applications, as it has an impact on the distribution of recharge to underlying bedrock aquifers.Modelling the geometry and distribution of the glacial deposits infilling the incised bedrock valleys with respect to their hydrogeologic properties allows the user to better understand groundwater migration pathways and the potential for increased groundwater recharge of bedrock aquifers.
Although the Paris Moraine (zone 2) in the subdivided regional model is dominated by till units, specific units of coarse-grained materials (outwash, kame, or ablation deposits) within the moraine tend to be thicker and more spatially constrained compared with that in the original regional model (Fig. 5).This is indeed more consistent with the subsurface stratigraphy of the moraine identified as heterogeneous and variable in regional studies based on water well records and geotechnical reports as well as more local studies based on both water well records and recovered sediment core (Blackport et al. 2009;Russell et al. 2013;Arnaud et al. 2017).Sediment composition in water well records of the Paris Moraine in the area surrounding Guelph, for example, suggests there is approximately 30% gravel, with relatively less sand (ϳ10%) and mud (<10%-15%) that make up the Paris Moraine subsurface stratigraphy (Russell et al. 2013).Cores recovered in the moraine together with cross section through the moraine in the Guelph area also have laterally restricted sequences of coarse-grained materials adjacent (within hundreds of metres) to boreholes dominated by diamict (Arnaud et al. 2017).The lateral continuity of units within zones 1 and 3 appear to have been more impacted by subdividing the model domain than those within the Paris Moraine.Generally, the geologic units within the subdivided regional model are more laterally discontinuous, and have greater variability in unit thicknesses (Fig. 5), which is also consistent with previous findings based on standard geological cross sections through the area (Blackport et al. 2009;McGill 2012;Arnaud et al. 2017), where subsurface heterogeneity, lateral discontinuity of geologic units, and (or) the relief on the underlying bedrock surface lead to rapid changes in unit thicknesses.For example, the basal diamict that underlies the outwash plain is not always preserved (McGill 2012), and the regional Catfish Creek Till thickens considerably from <5 m to over 25 m as it infills a topographic low in the underlying bedrock surface (Steelman et al. 2017).

Comparing range values
It is typically assumed that a longer range is better because it indicates that the data values were similar to one another over greater distances.However, in areas where the data are more variable over shorter distances, a variogram with a shorter range of correlation would likely produce a more realistic model result.The variogram ranges for the original regional model varied between 1818.9 m (outwash base grid) and 2496 m (ablation base grid), with an average range of 2210.8 m for the entire model (Table 2).The ranges for the subdivided models tended to be shorter than those of the original regional model (Tables 2-5).The range values in the zone 1 submodel varied between 1197.1 m for the outwash base grid and 1712.9 m for the undifferentiated till top grid, with an average range of 1585.9 m for the entire submodel (Table 3).The zone 2 submodel ranges varied between 1599.3 m for the Wentworth Till base grid and 2182.2 m for the undifferentiated till base grid, with an average range of 1973.8 m for the entire model (Table 4).The range values for the stratigraphic units within the zone 3 submodel varied between 1284.9 m for the Wentworth Till base grid and 1809.9 m for the kame base grid, with an average range of 1531.4 m for all of the stratigraphic units within this submodel (Table 5).
Glacial ice margin advances and retreats over the study area have led to a complex interfingering of glacial deposits of variable texture.Regionally, the subsurface stratigraphy has been characterized as consisting of alternating diamict and stratified sediment such as sand, gravel, and mud; the former represents ice advances, whereas the latter records ice margin retreat and glaciofluvial and glaciolacustrine processes in ice proximal to ice-marginal settings.However, local studies of the subsurface stratigraphy to date have demonstrated that the subsurface stratigraphy is more variable with laterally discontinuous stratigraphic units (Blackport et al. 2009;Russell et al. 2013;Arnaud et al. 2017).This heterogeneity is attributed to the dynamic interplay of glacial processes over time, namely glaciofluvial and glacial erosion and variable preservation of previously deposited glacial sediments, as well as the spatially and temporally varied distribution of glaciofluvial networks and glaciolacustrine ponding that would have occurred in different places over time as the ice margin melted and glaciofluvial and ice-marginal processes modified the preexisting deposits.
Considering the highly variable nature of the deposits described by various geological studies of the Paris Moraine (Karrow 1968;Barnett 1992;Sadura et al. 2006;Blackport et al. 2009;Russell et al. 2009;Arnaud et al. 2017), variograms with a shorter range are probably geologically more realistic, characterizing the true variability within the subsurface.By subdividing the study area prior to interpolation, the average range lengths for the stratigraphic units were decreased between 11% and 31% compared with the original regional model where all data were lumped rather than subdivided (Tables 2-5).A lower range value For personal use only.
does not suggest that one model is necessarily better than another, but rather that the variability has been constrained to the local submodel areas, and ensuring the stratigraphic variability characteristic to each submodel zone was preserved.

Comparing nugget values
The variogram nugget represents the amount of very short range (zero distance) variability that cannot be accounted for in the variogram model and is associated with unexplained errors in the data or sampling procedures (Krajewski and Gibbs 1996).The nugget is often the result of sampling errors or geological variability that is occurring on a smaller scale than the data point spacing and is measured in units squared (Davis 2002).Although it would be ideal to have a variogram without any nugget effect, this is typically uncommon in geological investigations (Krajewski and Gibbs 1996;Davis 2002).The nugget for the original regional model ranged from 0 to 25.3 m 2 for the seven stratigraphic units, with an average nugget of 14.9 m 2 for all units (Table 2).The zone 1 submodel had nugget values that ranged from 0 to 17.1 m 2 , with an average nugget of 11.5 m 2 (Table 3).The zone 2 submodel had nugget values that ranged from 0 to 15.6 m 2 , and an average nugget of 7.5 m 2 (Table 4).The zone 3 submodel produced variograms with nuggets ranging from 0 to 9.9 m 2 , and an average nugget of 4.7 m 2 (Table 5).Subdividing the regional model resulted in lower nugget values for each of the submodels, as average nugget values were reduced by 54%, 50%, and 68% for the zone 1, 2, and 3 submodels, respectively.These results suggest that

Comparing sill values
The sill represents the amount of variability observed between the data points at the variogram range.Theoretically, if the distance between the sample points is zero, the value at each data point will be compared with itself; therefore, the variance would also be zero.As the distance between points get slightly larger, so too would the difference in values.At some distance, the points are so far apart that they are no longer spatially related.At this point, the variance (gamma), which is measured in units squared, between successive data points no longer increases and it is referred to as the sill.If the sill is reached at a low gamma value, this indicates that there tends to be less variability between data point values.In contrast, if the sill is reached at a high gamma value, then there is generally more variability between data point values.Statistically, it is better to have a variogram with a low sill value, which indicates that there is lower variability between the data points even up to the distance at which the points are no longer statistically similar to one another (Krajewski and Gibbs 1996).However, geologically, a low or high sill value allows the level of variability within the stratigraphic dataset to be established, which in turn provides information on the subsurface variability within each of the submodel zones.The original regional model variogram produced sills ranging from 36.4 to 144.2 m 2 , with an average sill of 103.5 m 2 (Table 2).The zone 1 submodel variograms produced sills ranging from 10.6 to 108.3 m 2 , and an average sill of 76.6 m 2 (Table 3).The zone 2 variogram sill values ranged from 20.6 to 98.8 m 2 , with an average sill of 66.7 m 2 (Table 4).The zone 3 submodel variogram sill values ranged from 19.2 to 83.1 m 2 , with an average sill value of 45.2 m 2 (Table 5).The submodel variograms had lower average sill values than the original regional model (103.5 m 2 ; Table 2).A lower sill value does not suggest that one model is better than another, but rather that there is less variability within the data used to model the variogram.The results of this study indicate that there is less variability between the data points within the submodel variograms than in the original regional model variogram representing the entire model domain.The lower sill values are consistent with our current understanding of the local geology: zone 1 and zone 3 are characterized by the drumlinized till plain and outwash plain, and have relatively consistent subsurface stratigraphy dominated by till and stratified sediments, respectively, whereas the slightly larger sill value observed for zone 2 is representative of the more heterogeneous deposits of the Paris Moraine (Arnaud et al. 2017).By subdividing the regional model into submodels of similar data variability, the respective variograms were able to better characterize the site variability and are more likely to produce a more realistic representation of the volumes for each geologic unit.

Correlation coefficients
The correlation coefficient is used to assess how well the variogram model matches the data, which provides some information on the predictive capabilities of the interpolation parameters.The correlation coefficient is a measure between -1 and +1, indicating that there is a perfect negative correlation (-1), no correlation (0), or a perfect positive correlation (+1) between the data points and the variogram model.A correlation coefficient close to 1 indicates that the variogram was able to accurately define the spatial relationship of the data points and that there is a good fit between the observed data points and the variogram model that will be used to predict the grid surfaces; both of which are extremely important to achieve an accurate interpolation result.The correlation coefficient is therefore a useful statistic to help analyze how well the variogram model is able to match the observed data, which provides some information on the predictive capabilities of the interpolation parameters.In this study, we are interpolating the elevation of specific stratigraphic units; therefore, the correlation coefficient was used to assess how well the variogram model was able to predict the elevation of the stratigraphic units at observed (borehole) locations.The original regional model produced correlation coefficients between 0.89 and 0.99, with an average of 0.96 for all grids (Table 2).The zone 1 submodel produced variogram models with better correlation coefficients ranging from 0.98 to 1.0, with an average of 0.98 (Table 3).The zone 2 submodel correlation coefficients ranged between 0.97 and 1.0, with an average of 0.99 (Table 4).The zone 3 submodel values ranged between 0.88 and 0.99, with an average of 0.95 (Table 5).The zone 1 and 2 submodels both performed better than the original regional model, whereas the zone 3 submodel performed equally as well as the original regional model.
It is apparent that there is higher variability amongst the data of the zone 3 submodel due to the lower average correlation coefficient.This may be due to the complex depositional conditions that would have existed in this area during ice retreat and that include slumping of sediments along the Paris Moraine back slope as the buttressing ice retreated from the area, as well as glaciofluvial reworking and ponding that would have occurred as the ice retreated to the position of the Galt Moraine to the SE.In addition, some of the lower correlation coefficients observed in zone 3 may Pagination not final (cite DOI) / Pagination provisoire (citer le DOI) be related to the difficulty of identifying those units at depth (e.g., Port Stanley and undifferentiated tills; Table 5).This may be due to several factors, including (i) the possibility that these lower tills were partially eroded and therefore more randomly distributed, (ii) the predominance of low-quality data points in zone 3, which makes it difficult to distinguish between the various tills at depth, and (iii) the relatively sparse coverage of data points within portions of zone 3. It is likely that the lower correlation coefficient observed within the zone 3 submodel is causing the lower correlation coefficient observed within the original regional model, as both the zone 1 and 2 submodels produced variograms with better (higher) correlation coefficients (0.98 and 0.99, respectively; Tables 3 and 4) than the original regional model (0.96; Table 2).Therefore, by subdividing the regional model into submodel areas, it was possible to limit the amount of influence between the regions of increased variability observed within the zone 3 submodel, and the less variable model areas of zones 1 and 2. Developing a geologic model that more realistically characterizes the stratigraphy and subsurface heterogeneity will ultimately provide a more robust framework for assigning hydrogeologic parameters in groundwater flow models used in regional resource management or contaminated site remediation applications (de Marsily et al. 2005).

Optimal variogram model
For the variogram to be used in the interpolation process, a variogram model must be fit to the data to provide a prediction of the unknown values at specific locations.The variogram model is essentially a generalized mathematical equation that provides the best fit to the observed data (Isaaks and Srivastava 1989;Gringarten and Deutsch 2001;Chiles and Delfiner 2012).There are a variety of variogram models (i.e., spherical, exponential, linear, Gaussian, nugget, hole-effect, and drift) that can be applied to data, depending on how quickly the data variability increases as the distance between data points increases.Many software packages will automatically assign a spherical model when fitting a variogram unless otherwise specified.For this study, variograms were fit using each type of variogram model to determine which variogram model produced the best fit with the observed data points, had a high correlation coefficient, and corresponded with our understanding of the spatial variability for each specific stratigraphic unit (Gringarten and Deutsch 2001).A complete list of variogram models and their implications to model output can be found in Isaaks and Srivastava (1989), Krajewski and Gibbs (1996), Gringarten andDeutsch (2001), andJohnston et al. (2001).
The exponential variogram model most often provided the best fit for the stratigraphic unit datasets in this study.Exponential variogram models typically provide the best fit for datasets where the variance increases quickly at short distances and then gradually tapers off for data points that are further apart (Gringarten and Deutsch 2001).This seems appropriate considering many stratigraphic units in this study area are quite variable at short distances; however, because these units are often dispersed throughout the study area, the data tend to fit best with models with positive correlations that persist below the sill for large horizontal distances, and often do not reach the expected sill variance.All of the grid surfaces for the original regional model were modelled using an exponential variogram either with or without a nugget (Table 2).Only three surfaces in the zone 1 and 2 submodels (outwash deposit top in zones 1 and 2 and base in zone 1) were modelled with a Gaussian variogram model either with or without a nugget (Tables 3, 4).Gaussian variogram models typically provide the best fit for datasets whose variance initially increases gradually, then begins to increase more rapidly as the distance between points increases, which is likely considering the outwash deposits have a patchy distribution at the surface (Fig. 3C; Krajewski and Gibbs 1996).The only other grid that used a variogram model other than exponential was the ablation base grid in the zone 3 submodel, which used a spherical model (Table 5).Spherical variogram models are characteristic of datasets where variance initially increases linearly with increasing distance up to a welldefined sill, and are also typically used for modelling properties with higher levels of short-range variability (Krajewski and Gibbs 1996).In general, there was very little change in variogram type between the original regional model and the individual submodels; therefore, subdividing the regional model domain for this study area was minimally impacted by the type of variogram selected.

Variation in the direction of longest variogram range
To determine the direction for which the variogram range is the longest and produces the greatest correlation over the greatest distance, variograms were evaluated in multiple directions to test the data for anisotropy.The direction of longest variogram range provides a sense of the direction in which sediment types are most consistent from one data point to the next.This parameter may change when the regional model is subdivided, if dominant processes change from one area of the regional model to the next.The direction of longest variogram range is also expected to be more closely aligned for the top and base of geological units in a submodel zone, as that zone will have been subjected to specific depositional conditions and be associated with specific geological features that are directionally more consistent, compared with the wide range of depositional conditions and anisotropy captured in a regional model.Lower correlation between the top and base of thicker or more extensive units is expected where the basal contact and the overlying sediment result from different processes (e.g., glaciolacustrine sediments blanketing an irregular topography associated with the underlying depositional regime).
The results from the anisotropic variograms revealed that the direction of longest range changed for some stratigraphic units when the regional model was subdivided into submodels based on data variability and geological context.For example, the top and bottom grid surfaces of the kame deposit for the original regional model produced the longest ranges at N135°and N141°, respectively.When the variograms were produced for the subdivided models, the top and bottom grid surfaces of the kame deposits for the zone 2 submodel produced the longest ranges at N7°and N9°, respectively, which represents a 130°shift (Tables 2, 4).Kames can have a wide variety of ranges considering that they are deposits that represent the internal drainage that was occurring on and within the ice (Benn and Evans 2010).The change in the direction of longest range between the original regional model and the zone 2 submodel highlights the variation of the kame deposits between the submodel zones while maintaining consistency within each submodel zone.
The consistency of the range direction did indeed increase for two of the three submodel regions when the regional model domain was subdivided.The direction of longest range within the original regional model varied between N3.9°and N169°(degrees from north; Fig. 6A; Table 2) for the seven stratigraphic units.The direction of greatest correlation range for the submodels was much more consistent within the zone 1 and 2 submodels.The major axis for stratigraphic units within zone 1 ranged from N121°t o N149°, which is a difference of only 28°, and coincides with the orientation of drumlins, the direction of ice flow, and meltwater runoff in this study area (Karrow 1968(Karrow , 1987;;Barnett 1992; Fig. 6B; Table 3).There was also much greater consistency between the top and base grids for each stratigraphic unit within zone 2 (average difference of only 12.6°; Fig. 6C; Table 4).Although the range direction for the zone 2 submodel stratigraphic units varied between N7°and N176°(a 169°difference), the rose diagram for the zone 2 submodel (Fig. 6C) shows that the direction of longest range for the majority of stratigraphic units were aligned in the direction of ice movement and meltwater runoff to the NW.In contrast, the kame deposits and uppermost outwash deposits of Pagination not final (cite DOI) / Pagination provisoire (citer le DOI) zone 2 were aligned nearly N-S (Fig. 6C; Table 4), which in some places is aligned with the front edge of the Paris Moraine (Fig. 1) and explains the 169°difference in range directions for the different stratigraphic units of zone 2. It is perhaps not surprising that kames are not consistent with the rest of the data considering their deposition is related to glacial hydrology controlled by localized ice configuration and drainage networks.In contrast, the zone 3 submodel showed considerable variation in the direction of the unit variogram long axis.The directions ranged from N8.6°t o N177°, which is a difference of 168.4°(Fig.6D; Table 5).This is consistent with the paleoenvironmental model for this area that suggests initial deposition and deformation by ice in the SE-NW direction, further deposition by debris flows as the ice retreated and support of the moraine sediment was removed, and lastly glaciofluvial and glaciolacustrine reworking of sediments along the back slope of the Paris Moraine in the NE-SW direction as water drained between the Paris Moraine and the adjacent Galt Moraine.The difference in direction between the top and base grid surfaces for each unit was also significantly lower, with an average difference of 54.4°(Fig.6D; Table 5).
Visually, the 3D stratigraphic models of the original regional model and the subdivided regional model appear quite different from one another (Figs.3-5).The most substantial difference is in the outwash deposit, which appears much less extensive in the subdivided regional model, and is limited to two broad swaths; one in zone 1 and one in the SE corner of the study area in zone 3 (Fig. 3).Zone 2 of the subdivided model is characterized by more extensive tracts of Wentworth Till, kames, and minor outwash Pagination not final (cite DOI) / Pagination provisoire (citer le DOI) (Fig. 3B), which are consistent with the distribution of Wentworth Till and gravel identified in the surficial geology map (Fig. 3C).The variation in the range direction between the top and bottom of the outwash stratigraphic unit in the original regional model was 40°(Fig.6A; Table 2).However, when the outwash unit was interpolated in the submodel sections, the variation between the top and bottom grids was a maximum of 20°(Tables 2-5).The reduced variation in range likely occurred because there was greater correlation between the direction of longest range between the top and bottom grids of the stratigraphic units, which resulted in a more localized and streamlined model representation of the outwash deposits (Fig. 3).In contrast, the deposits within zone 3 had much more directional variability than those of zone 2 or 1 (Fig. 6D).It is likely that the greater directional variability of the stratigraphic units within zone 3 significantly affected the original regional model variograms, resulting in more directional variation and greater inconsistencies between top and bottom of each unit, therefore biasing the form and geometry of the outwash unit.By subdividing the regional model, it was possible to produce variograms that better characterize the individual stratigraphic units within each zone, thus revealing the different degrees of variability associated with each zone, reflective of depositional environments and spatial data.
In geological terms, the subdivided model variograms were able to capture the variability in flow direction associated with each zone and thus better infer the resulting distribution of the outwash stratigraphic unit in those different meltwater drainage networks.In zone 1, the direction of longest variogram range is consistent with meltwater flow directions away from the ice margin and the Paris Moraine and influenced by the NW-SE drumlin orientations of the preexisting till plain.In zone 2, the direction of longest variogram range for the top and base of the outwash deposit varies by 20°, consistent with the likely more variable drainage patterns influenced by the chaotic hummocky topography of the moraine.Lastly, in zone 3, the direction of longest variogram shifts to a NNW-SSE direction, which is perhaps capturing some of the predominant NW-SE drainage as well as some deflection of that flow towards the south as meltwater was routed between the Paris and adjacent Galt moraines.

RMSE results
Root mean square error results have been used in many studies as an effective method of assessing the accuracy (i.e., how closely the model matches the observed points) of interpolated results (Zimmerman et al. 1999;Davis 2002;Dille et al. 2003;Mueller et al. 2004).The RMSE results were used in this study to compare the model accuracy of the regional model interpolated using all the available data (original regional model) and the regional model produced by subdividing and interpolating submodels based on data complexity and (or) variability (subdivided regional model; Fig. 3).The subdivided regional model produced a lower RMSE value (0.57) than the original regional model (0.86).Therefore, subdividing the regional model domain into submodels produced a more accurate spatial interpolation of the regional data.
Ideally, only data from high-quality and reliable data sources would have been used.However, as is often the case in regional geological models, it was necessary to include variable quality data from regional databases owing to the size of the regional model domain, and the fact that the available moderate-quality data were often clustered and did not penetrate very deep in the subsurface.It is likely that as additional high-quality data become available for the model domain, the performance of the subdivided model variograms will improve further, thus having an even greater impact on the 3D stratigraphic model accuracy and an improved basis for informing hydrogeologic parameter variability.Stochastic modelling approaches could also be applied within submodel zones to visualize and assess multiple model realizations and capture the variability in the model predictions (Kearsey et al. 2015).The next step is to integrate this improved representation of the geologic spatial variability with hydrologic parameters to assess the effects on flow system conditions where water balance and hydraulic head distributions from field and model outputs can be evaluated.This improved delineation of the subsurface deposits within the Paris Moraine will in turn be useful for the management of groundwater resources as source water protection policies are implemented in coming years.

Conclusions
The objective of this study was to determine if subsurface geological model accuracy could be improved by splitting the regional model into multiple submodels based on the degree of variability observed between surrounding data points and our knowledge of the geological-depositional framework from previous studies.The hypothesis is that by separating the regional model domain into submodels based on the data variability and our understanding of the geological context, a variogram can be produced for each submodel region that is more representative of variability within each area, thus resulting in a more geostatistically appropriate and geologically realistic (cf.Karrow 1968) model interpolation for the entire study area that can be used to better inform aggregate and groundwater resource characterization studies (Weissmann et al. 1999;de Marsily et al. 2005).
The results showed that by interpolating the regional stratigraphic units within submodels based on the spatial variability within the data, the kriging algorithm was able to produce a more representative variogram of the data.The zone 1 and 2 submodels produced variograms with lower ranges, more consistent major axis directions, higher correlation coefficients, and lower sill and nugget values.This is to be expected considering the stratigraphic units in zones 1 and 2 are typically elongate and laterally continuous deposits in the SE-NW direction for the outwash and drumlinized till plain deposits (zone 1), recording the dominant meltwater and ice flow directions, and in the NE-SW direction along the Paris Moraine (zone 2).Thus, the variogram models for these submodels were more constrained and able to better characterize the observation data than when the data for the entire model domain were considered as a whole.The zone 3 submodels produced variograms with the shortest ranges, lots of variability in the major axis direction, and produced lower correlation coefficients.The zone 3 submodel data were evidently more difficult for the variogram to characterize, suggesting that this area contains the most variable and complex stratigraphy within the model domain.This variability is attributed to sediment reworking by gravity and glaciofluvial processes, and drainage water ponding between the topographic highs of the Paris Moraine and nearby Galt Moraine.Zone 3 likely requires additional high-resolution field data to help characterize this complexity or further subdivision to improve variogram parameters and model outputs.
The results of this study reveal that when all of the data were modelled together as one large regional model, the resulting variograms produced blended estimates of the variogram parameters.However, when the zone 3 submodel data were modelled independently of the zone 1 and 2 submodels, the variograms were better able to capture the subsurface characteristics of each zone.The increased variability within zone 3 attributed to the complex sediment reworking and variable glaciofluvial, glaciolacustrine, and ice-marginal processes associated with ice margin retreat to the adjacent Galt Moraine was confined to zone 3 and did not impact the submodels of the neighbouring zones.Conversely, the data characteristics of the more laterally extensive and directionally consistent stratigraphic units within the outwash and till plain of zone 1 as well as the elongate form and geometry of the stratigraphic units within the Paris Moraine (zone 2) submodels were confined to those areas and did not reduce the variability observed within the zone 3 submodel area.
Pagination not final (cite DOI) / Pagination provisoire (citer le DOI) The proposed submodel approach produces a more geologically realistic model of the study area when compared with the regional model approach as evidenced by comparing the model results against the regional surficial geology map (Ontario Geological Survey 2003) and previous subsurface studies in this area (Karrow 1968(Karrow , 1987;;Blackport et al. 2009;Arnaud et al. 2017).In the submodel outputs, various stratigraphic units are laterally more constrained, and their bounding surfaces, including those of the underlying bedrock, exhibit more significant relief.This is consistent with subsurface observations in the study area where tills are laterally discontinuous, stratified sediments are found in close proximity to diamict within the moraine, and the underlying bedrock surface is dissected by various well-defined buried bedrock valleys.The submodel results are also more consistent with what is generally known of sediment landform associations in terrestrial glaciated landsystems, namely the relative lateral continuity of facies in till and outwash plains and the relative variability of facies found in ice-marginal or ice-contact sediments (Stephenson et al. 1988;Benn and Evans 2010).The integration of submodels provides a means to highlight the heterogeneity across the regional domain and avoids blending areas of differing variability together, resulting in a more realistic characterization of the nature, geometry, and extent of geological and hydrogeological units within this complex depositional environment.This is turn can improve the geological framework used in various groundwater flow modelling efforts, both in site-specific contaminant transport investigations or in regional susceptibility or groundwater supply studies carried out in previously glaciated regions like southern Ontario.

Fig. 1 .
Fig. 1.Map of the regional model domain, data point locations, and locations of cross sections (A-A 1 , B-B 1 , and C-C 1 ).There are 1245 low-quality data points (red circles) from the water well database, 263 moderate-quality data points (yellow squares) from more reliable sources such as cored boreholes and geotechnical reports, and six high-quality data points (green triangles) from related research of the co-authors.The three submodel regions are delineated with white lines and identified as zone 1, zone 2, and zone 3. Note that zone 2 broadly coincides with the Paris Moraine.Coordinates are in Universal Transverse Mercator (UTM) zone 17N; datum: North American Datum 1983 (NAD83).[Colour online.]

Fig. 2 .
Fig. 2. Three-dimensional borehole representation of the data distribution and regional stratigraphy.The vertical exaggeration of this model is 20×.[Colour online.]

Fig. 4 .Fig. 5 .
Fig. 4. Fence diagram through (A) the original regional model and (B) the subdivided regional model.Both fence diagrams have a vertical exaggeration of 20×.[Colour online.]

Table 2 .
Variogram parameters for the original regional model.

Table 3 .
Variogram parameters for the zone 1 submodel.

Table 4 .
Variogram parameters for the zone 2 submodel.

Table 5 .
Variogram parameters for the zone 3 submodel.