Accurate inference of the interfaces between geological units showing different hydraulic properties is a key step for a reliable hydrogeological characterization of regional aquifers. In this study, we developed a workflow that combines multiple geological and geophysical data sets having different intrinsic resolution to map a stratigraphic interface of the regional aquifer located in Montérégie, Quebec, Canada. One of the principal goals was to optimally assimilate all the data at all stages in the workflow. Firstly, the experimental variogram showed two structures of different ranges: one coming from highly sampled geophysical data and the other from conventional borehole geological markers. Secondly, a secondary variable is constructed with all secondary data (e.g., geological interpretations of low-resolution electromagnetic surveys), each having its own accuracy and resolution. To account for the variable secondary data, different reliability indexes were assigned as weights in a discrete smooth interpolation (DSI). Thirdly, a classic kriging with an external drift (KED) operator was used to interpolate the more reliable well data on the entire region. The approach was tested on the estimation of the bedrock interface elevation in a regional hydrogeological characterization study. The resulting map shows bedrock elevations coherent with geological structure of the region, representing main features such as outcrops and valleys. A cross section is presented to illustrate the philosophy behind the tools employed to achieve the estimation process. It also shows an example of visual quality control undertaken to validate the workflow.