A Python Application for Visualizing the 3D Stratigraphic Architecture of the Onshore Llobregat River Delta in NE Spain

Abstract

This paper introduces a Python application for visualizing the 3D stratigraphic architecture of porous sedimentary media. The application uses the parameter granulometry deduced from borehole lithological records to create interactive 3D HTML models of essential stratigraphic elements. On the basis of the high density of boreholes and the subsequent geological knowledge gained during the last six decades, the Quaternary onshore Llobregat River Delta (LRD) in northeastern Spain was selected to show the application. The public granulometry dataset produced by the Water Authority of Catalonia from 433 boreholes in this strategic coastal groundwater body was clustered into the clay–silt, coarse sand, and gravel classes. Three interactive 3D HTML models were created. The first shows the location of the boreholes granulometry. The second includes the main gravel and coarse sand sedimentary bodies (lithosomes) associated with the identified three stratigraphic intervals, called lower (>50 m b.s.l.) in the distal LRD sector, middle (20–50 m b.s.l.) in the central LRD, and upper (<20 m b.s.l.) spread over the entire LRD. The third deals with the basement (Pliocene and older rocks) top surface, which shows an overall steeped shape deepening toward the marine platform and local horsts, probably due to faulting. The modeled stratigraphic elements match well with the sedimentary structures reported in recent scientific publications. This proves the good performance of this incipient Python application for visualizing the 3D stratigraphic architecture, which is a crucial stage for groundwater management and governance.

1. Introduction

The Llobregat River Delta (LRD) is a densely populated coastal plain of 98 km², forming the southwestern sector of the metropolitan area of Barcelona city in the Catalonia region in northeastern Spain (Figure 1). This area also includes other minor-order cities such as El Prat de Llobregat, L’Hospitalet de Llobregat, Cornellà, Sant Boi, Viladecans, Gavà and Castelldefels. The permanent surface-water provision and the abundance of groundwater have favored the development of coastal aquatic ecosystems and human settlements since ancient times. During the XIX and XX centuries, a large industrial activity developed in the area. Since the second half of the XX century, the high groundwater exploitation rates to supply to the increasing population and industrial activity produced negative consequences on groundwater quantity and quality, including seawater intrusion into aquifers and high levels of pollution [1]. Other modern development milestones occurred in Barcelona city and its metropolitan area, such as the Olympic Games in 1992 and the Llobregat Delta Infrastructure and Environment Plan (Pla d’Infraestructures i Medi Ambient del Delta del Llobregat, PDL) started in 1994 [2] have modified the LRD land use and geomorphology. The PDL included large civil infrastructures with variable underground development such as the extension of the International Airport and Harbor of Barcelona, a dense network of new highways, different underground lines, conventional and high-speed railways, wastewater new treatment plants, desalination plants, diversion of the Llobregat River final section, and many other civil works.

The PDL was a disruptive milestone in what concerned the generation of new, testable geological and hydrogeological information. The conventional geotechnical exploration data for designing civil work foundations and the drilling of dozens of boreholes to monitor the impacts of civil works on groundwater provided huge, detailed geological and hydrogeological data. In 2004, the Water Authority of Catalonia (Agència Catalana de l’Aigua, ACA) created the Technical Unit of the Llobregat Aquifers (Mesa Tècnica dels Aqüífers del Llobregat, METALL). This Hydrogeological Research Unit was devoted to compiling and homogenizing the former and new (in production) geological and hydrogeological information in order to prepare groundwater flow numerical models aimed at assessing the cumulative impact of the large civil works on the groundwater bodies [3]. Nowadays, the ACA manages this public geological and hydrogeological database, which is available on request.

From a geological point of view, the LRD is regionally fed by the Llobregat River and its tributaries arriving from the Pre-Pyrenean Range and locally from the Llobregat River lower valley reliefs (Garraf and Collserola massifs) belonging to the Catalonian Coastal Range [4,5]. This range is a NE–SW-oriented mountain chain that gives pass downward to the Mediterranean coast (Figure 1). In the LRD, the geological studies started at the end of the XIX century by Almera [8]. In the XX century, several studies [4,5,6,9,10] allowed proposing geological maps and 2D cross-sections aimed to support hydrogeological studies. Sedimentological studies performed in the 1970s and 1980s [11,12,13] allowed clarifying the geology of the LRD. In the 1980s, the prodeltaic bodies of the emerged delta, dated as Holocene, were studied in detail [14] and the geological characterization of the continental shelf with support of 2D marine seismic reflection took place [12,13,15]. These studies, with a strong sedimentological component, allowed the sequential division of the LRD and the arrangement of the Quaternary materials. In the 1990s, this huge geological background was combined with the former geological information compiled from dozens of boreholes to make modern groundwater evaluations [16]. Coinciding with the PDL development at the beginning of the XXI century, new geological data allowed fine research [7,17,18,19,20] aimed to detail the 2D stratigraphic architecture of the LRD, thus giving a response to the Pliocene–Quaternary boundary, the confident definition of those detritic levels officially cataloged as productive, high-yielding aquifers and the detailed identification of those interconnected stratigraphic structures (of sedimentary and tectonic origin) through which seawater intrusion into aquifers and mobilization of pollutants take place.

In recent times, numerical modeling has progressively been employed to characterize the geometry of geological structures and sedimentary bodies (lithosomes). Numerical tools involve a variety of deterministic and stochastic models, as well as geostatistical and geospatial techniques. The choosing criteria depend on the purpose of the research, input data typology and spatial coverage, and capability for calibrating and validating predictions. Regarding the interpolation tools, this field of research allows the 3D mapping of certain physical parameters controlled by sedimentary processes that determine the hydraulic behavior of the geological reservoir. The 3D mapping allows better visualizations than the classic 2D serial cross-sections [21,22,23,24]. A variety of software for 3D visualization based on different interpolation algorithms is available. Commercial tools such as MOVE (Petroleum Experts Ltd., Edinburgh, UK), 3D Geomodeller (Intrepid Geophysics, Brighton, Australia), Autocad Civil (Autodesk, Inc., San Rafael, CA, USA), Gocad (Emerson Paradigm Roxar), ArcGis, PETREL (Geology and Modeling from Schlumberger), VOXI (Earth Modeling from Geosoft) and Geoscene 3D (I-GIS) are available for 3D geological visualizations. There are other open-source applications such as Gempy or OSGeo. The Gempy applications have proven to be accurate enough and have technical support for users, but they are expensive. The advantage of OSGeo applications is the zero cost and adaptability (modifying or extending the sources), but the absence of technical support for users and sometimes low reliability are its negative parts. There are commercial and open-source libraries devoted to geographic information systems (GIS) and mapping. For instance, there are Python libraries [25,26,27] and posts listing libraries [28,29] for different geological interests, including structural geology, sedimentology, mining, and hydrogeology. There are also scientific documents where computer tools are developed or applied to different fields of geology [30,31,32,33]. In addition, some social media channels have also posted data analytics and machine learning educational applications focused on geology [34].

The experience of the mathematician researchers involved in this paper with Python for data handling and 3D mapping was an advantage for adapting this knowledge to geological concepts and concerns. This adaptation required less effort than learning specific tools developed for geologists. So, this paper introduces a Python application for visualizing the 3D stratigraphic architecture. The application interpolates the parameter granulometry deduced from borehole lithological records to create interactive 3D HTML models of essential stratigraphic elements. For this, the application uses different Python packages, including basic interpolation and shape algorithms. On the basis of the high density of boreholes and the subsequent geological knowledge gained during the last six decades, the Quaternary onshore LRD was selected to show the application. The public geological and hydrogeological database managed by the ACA was consulted to compile the information used in this paper, i.e., the granulometry dataset generated by METALL.

Different research disciplines in Earth and Environment Sciences require 3D visualization of available geological information. As described above, the LRD is a strategic coastal aquifer whose sustainability is threatened by strong human pressure. In this case, the Python application may assist to drill proper pumping wells and optimize groundwater monitoring networks, as well as evaluating how the large civil works may impact the groundwater resources. The spatial distribution of granulometry will also be of assistance to design compensatory measures for aquifer protection and recoveries such as artificial aquifer recharge and positive hydraulic barriers to reduce the advance of seawater intrusion and mobilization of pollutants.

Appendix A includes the four Jupyter notebooks describing the data handling and 3D mapping. An operative version of the introduced Python code can be downloaded from the GitHub repository described in Supplementary Materials. The results of this Python application are HTML files that do not require any additional tools apart from a web browser to view different perspectives, hide elements, enlarge or focus on specific areas or elements and take snapshots of a particular view. These actions can all be performed using just a mouse.

2. Geological and Hydrological Setting

The LRD is located at the foot of the Catalonian Coastal Range (Figure 1). In the area, this mountain range includes rocks of Paleozoic (granites and slates) and Mesozoic (Triassic conglomerates, sandstones, and pelites; Jurassic dolostones and limestones; Cretaceous marly limestones). The LRD is NE bounded by the Montjuïc Mount (Figure 1) made of Miocene calcarenites and marly limestones [17]. This area represents a Neogene rifted margin associated with the opening of the Valencia Trough and is affected by several fault families (Figure 1) probably active and mainly oriented NE–SW (Morrot and Tibidabo fault families) and NW–SE (Llobregat fault family), that conditioned the main reliefs and the location of the Llobregat River outlet [7].

With respect to the LRD formations (Figure 2), Pliocene and older rocks are considered the basement, which is separated from the Quaternary formations by an important unconformity surface [18]. Pliocene substratum is made of estuarine marls, silts, and clays [12,13,17]. The Quaternary record was divided into two depositional sequences Pleistocene and Holocene in age [7,12,13,17,18,20]. The terms upper detrital complex and lower detrital complex have also been used [7] for the same depositional sequences (Figure 2). According to geophysical studies in the offshore delta in the marine platform [12,13], the lower detrital complex can be divided into three parasequences [7]. In general terms, the lower detrital complex is made of conglomerate bodies (locally with sand) with intercalated silt- and clay-rich intervals [12,13,35]. The upper detrital complex, from bottom to top, is made of a sand layer, a silt bed, gravel (locally with sand), and upper silt and clay cover forming the current alluvial plain and the associated coastal wetlands and marshes [11,35].

From a hydrological point of view, the Llobregat River basin spans from the Pyrenees (the mountain range located between Spain and France; Figure 1) to the Mediterranean Sea with a total length of 170 km and covering around 4950 km² [20]. The Llobregat River shows an irregular streamflow regime in the 250–1500 Mm³ range determined by a Mediterranean-like climate regime [11]. The deltaic plain is also modeled by several ephemeral streams coming from the Garraf and Collserola massifs [17]. Apart from this hydrological control, the LRD coastal fringe has also been modified by the regional littoral drift towards the SW. In this area, littoral currents of about 30 cm s^–1 redistribute the sediments to the SW also bringing sediments from northern rivers and streams [36,37]. To a lesser extent, the waves (low energy) and the tide (a few centimeters) contribute to the sediment’s redistribution in the coastal areas of the LRD [36,37].

Figure 2. Geological sketch cross-sections of the LRD located in Figure 1 ((A-A’) is SW–NE oriented; (B-B’) is NW–SE oriented), modified and updated from MOP [35], Marqués [11], Simó et al. [20], Gámez et al. [7].

Figure 2. Geological sketch cross-sections of the LRD located in Figure 1 ((A-A’) is SW–NE oriented; (B-B’) is NW–SE oriented), modified and updated from MOP [35], Marqués [11], Simó et al. [20], Gámez et al. [7].

3. Materials and Methods

3.1. Python Programing Language

The object-oriented, high-level Python programming language [38] was used for granulometry data treatment and displaying the 3D HTML models of essential stratigraphic elements. Python is an open-source language used in many environmental topics, including geological ones. Many packages or modules use this programming language for a wide variety of problems. The Python packages used here were (i) NumPy [39] for data computing, (ii) Pandas [40] for data analysis and processing, (iii) Plotly [41] as a plotting library, (iv) UTM [42] to deal with UTM coordinates and borehole prospecting depths, and (v) Scipy [43] for interpolation and render algorithms.

3.2. Data Compilation and Pretreatment

The public geological and hydrogeological database managed by ACA was consulted to compile the lithological records of 433 boreholes in the onshore LRD, which can be also downloaded from Alcalá et al. [44], and the granulometry dataset generated by METALL from both Lab tests values and proxy values after visual recognitions of the predominant lithologies, which is available on request. The geotechnical borehole SM P-4 (located in Figure 1) with granulometry data is included in Figure 3 as an example. The dataset consisted of georeferenced XLS (Excel) files with meter-by-meter granulometry values associated with the boreholes lithotypes. The boreholes’ location (coordinates x and y) and prospecting depth (coordinates z) lead to a georeferenced array of granulometry data associated with the prospected lithologies over space and depth (Figure 4).

Figure 3. A selected geotechnical borehole (SM P-4, location in Figure 1), showing the granulometry tests (red arrows) performed through the Spanish standards and their results for gravel (cyan), coarse sand (yellow), and silt–clay (gray) granulometry classes. This information forms part of the granulometry dataset generated by the METALL used in this paper.

Figure 3. A selected geotechnical borehole (SM P-4, location in Figure 1), showing the granulometry tests (red arrows) performed through the Spanish standards and their results for gravel (cyan), coarse sand (yellow), and silt–clay (gray) granulometry classes. This information forms part of the granulometry dataset generated by the METALL used in this paper.

The georeferenced granulometry dataset and the original boreholes lithologies [44] were additionally compared to detect possible outliers. After checking, the georeferenced granulometry data were clustered into three main granulometry classes: clay–silt (<1 mm), coarse sand (1–5 mm), and gravel (>5 mm). Although expandable because Pandas can handle XLS files directly, data conversion to CSV format for further processing with Pandas has been preferred. The CSV files are easier to handle using a plain text editor only, so the users can use Pandas or just a plain text editor to make data edits needed.

3.3. The 3D Mapping of the Boreholes Granulometry Classes

Python has different visualization libraries, and Plotly [40] is one of them. Plotly was chosen because it is a web-based tool kit that can be accessed from a Python Notebook. The Plotly function Scatter3d was used for mapping boreholes. The use of this function requires additional data processing.

First, the function coordinates (data, positions) were defined. This function (i) lists the X, Y, and Z UTM coordinates extracted from ‘data’ by looking at the data indicated at ‘positions’, and (ii) models the data in legible format by the Scatter3d function. However, only the boreholes inside the LRD contour are needed. The function near (xyz,polyg,dis) is defined to select these boreholes. This function uses the Python Geometry function ‘distance’ to select coordinates within a distance less than ‘dis’ from the polygon ‘polyg’ in the ‘xyz’ list. Once these data are obtained in the required format, the Plotly function Scatter3d is used to create the figure. Some parameters such as the symbol, size, and opacity of the borehole, thickness and type of the LRD contour, and figure limits must be chosen. The aspect ratio determines the proportion of the Z axis regarding the X and Y axes, thus allowing the exaggeration of the Z scale for better displaying. This ratio must also be defined. This working flow is synthesized in Figure 5. Appendix A and the Jupyter notebook boreholes_3D_map.ipynb hosted in the GitHub repository described in Supplementary Materials include the complete code.

3.4. The 3D Mapping of the Stratigraphic Architecture

Due to the hydrogeological interest of the LRD, we decided to focus the study on the gravel and coarse sand lithosomes forming the local aquifers (Figure 6). When checking the borehole lithosomes (3D_Boreholes_LRD.html), gravel and coarse sand lithologies appeared concentrated preferably in certain zones of the top, middle and low levels. Then, a list of points (q_up, q_mid, and q_low for gravel, and p_up, p_mid, and p_low for coarse sand) with a representative point in each zone with a predominant granulometry was defined. A mathematical strategy (nucleation-like) is further used to create (and shape) the lithosomes associated to a given granulometry class (quasi-equal granulometry lithosomes). The mathematical strategy involves a function named grouping that depends on two auxiliary functions named within and within2. The function within (p,list,n) depends on three parameters: a point ‘p’, a list of points ‘list’, and a number ‘n’. This function is defined in plain Python (does not need to use NumPy functions) and just selects the elements of the variable list located within a distance n of the point p. The other within2 (list1,list2,n) function, also defined in plain Python, applies the function within (p,list2,n) for each p in list1 and puts the selected elements in the return list. Finally, the grouping (list1,list2,n) function recursively uses the within2 function to return a list with all elements in the variable list2 located within distance n of any element that is itself within distance n of each given point in the variable list1. This third function is also defined in plain Python.

Clearly, the efficiency of this function grouping is improvable, since its recursive nature implies repeating the same calculation many times, but its definition is very simple, and it does the job. Nevertheless, the lack of information between boreholes forces us to manually infer the lithological type in those blank spaces. To solve this gap, we add two lists (qq for gravel and pp for coarse sand) to the original data at each level. This makes the grouping strategy work properly.

Once the granulometry data of each lithosome were grouped by using the above nucleation strategy, the elements forming the previously defined groups of each lithosome were computed. For this, the Convex Hull algorithm developed by the SciPy community [44] was implemented. The 3D convex hull of a georeferenced dataset is the smallest polyhedron that wraps up all of them. The convex hull function must be applied to each obtained lithosome. The Delaunay triangulation algorithm [45] is used to render the shape of the lithosomes.

Finally, the ‘Mesh3d’ Plotly function was implemented to draw the volume figure, named 3D_Lithosomes_LRD.html. The output is an interactive 3D HTML model that (i) can be opened with any browser and (ii) allows observing different views of the mapped lithosomes, zooming, rotating, and moving around, as well as hiding any of them by clicking in the proper element in the right-side legend. This working flow is synthetized in Figure 6. Appendix A and the Jupyter notebook boreholes_3D_map.ipynb hosted in the GitHub repository described in Supplementary Materials include the complete code.

3.5. The 3D Mapping of the Basement Top Surface

A basic inverse distance weighted (IDW) interpolation algorithm was implemented to map the basement top surface (Figure 7). This IDW spatial function was obtained from the GEODOSE block, which can be found in the 3D terrain modelling package on the Python website [46]. To run this spatial function, a maximum data searching over a maximum searching distance (in both x, y, and z coordinates) was imposed, taking the boundaries imposed by the onshore LRD contour (x and y coordinates) and basement depth (coordinate z is the basement first-prospecting depth) into account. The first step was selection of those boreholes reaching the basement within the onshore LRD contour. This information was gathered from the CSV files to predict (by using the IDW spatial function) the continuous basement top surface shape.

The mapping grid has to be chosen in such a way that the computing time will be reasonable, and the result obtained was geologically sound. For this, it is important to choose a proper iteration radius (block radius) and the p-value ‘p’ for the weight that the IDW spatial function uses. Using the selected mapping grid, the Plotly function surface is implemented to draw the basement top surface.

The output is an interactive 3D HTML model that can also be opened with any web browser and allows observing different views and details of the mapped basement top surface. This working flow is synthetized in Figure 7. Appendix A and the Jupyter notebook boreholes_3D_map.ipynb hosted in the GitHub repository described in Supplementary Materials include the complete code.

4. Results

4.1. The 3D Mapping of the Boreholes Granulometry Classes

Figure 8 shows the 3D HTML model created to display the 3D distribution of the granulometry classes along the Z axis in each of the 433 boreholes complied in the onshore LRD. An interactive 3D HTML version of this model (3D_Boreholes_LRD.html) is provided as supplementary material. The plotting adopted a 1:1:50 (x = 2, y = 2, and z = 0.5) aspect ratio because of the thinness of the LRD Quaternary record regarding the x and y distances. This 3D view allows identifying three intervals where coarse detritic rocks (gravel and coarse sand) are concentrated and separated by other intermediate fine detritic intervals (silt and clay) as (i) >50 m b.s.l. lower interval (b.s.l. means below sea level), (ii) 20–50 m b.s.l. middle interval, and (iii) <20 m b.s.l. upper interval. Moreover, this 3D view also allows defining several bodies grouping the main gravel and coarse sand lithosomes. These stratigraphic intervals including the lithosomes were proposed by standard stratigraphic correlation procedures.

4.2. The 3D Stratigraphic Architecture of the Quaternary Coarse Detritic Lithosomes

Figure 9 shows the 3D HTML model created to display the 3D distribution of the gravel (G) and coarse sand (S) lithosomes within the onshore LRD lower (l), middle (m), and upper (u) intervals. An interactive 3D HTML version of this model (3D_Lithosomes_LRD.html) is provided as supplementary material. The plotting adopted a 1:1:50 (x = 2, y = 2, and z = 0.5) aspect ratio to show three different zenithal, SE–NW frontal, and NE–SW lateral frontal views of the lithosomes, which are labeled as Gl-1 and Gl-2 for the gravel lithosomes of the lower interval, Sl-1 for the coarse sand lithesome of the lower interval, Gm-1 to Gm-5 for the gravel lithosomes of the middle interval, Sm-1 and Sm-2 for the coarse sand lithosomes of the middle interval, Gu-1 to Gu-5 for the gravel lithosomes of the upper interval, and Su-1 to Su-5 for the coarse sand lithosomes of the upper interval. The dark gray lines are rendering effects of the implemented Delaunay triangulation algorithm to improve the shape of the lithosomes volumes and highlight the faces of the convex hulls.

4.3. The 3D Mapping of the Basement Top Surface

The basement top surface is an important stratigraphic element that determines the accommodation space for the Quaternary sedimentation. Within the onshore LRD contour, a total of 87 boreholes reaching the basement were gathered from the CSV files to predict (by using the IDW spatial function) the continuous 3D basement top surface. Taking the 98 km² onshore surface of the LRD into account, after several trials, the optimal regular mapping grid was 100 m × 100 m (n = 100), the block radius was 5 (r = 5), and the p-value was 2 (p = 2). Figure 10 shows the 3D HTML model created to display the 3D mapping of the basement top surface. An interactive 3D HTML version of this model (3D_Basement_LRD.html) is provided as supplementary material. The plotting adopted a 1:1:50 (x = 2, y = 2 and z = 0.5) aspect ratio to show three different zenithal, SE–NW frontal and NE–SW lateral frontal views. This model also allows different views, zooming, rotating, and moving around, as well as hiding elements to focus on details.

5. Discussion

From the above-introduced methodology and the scope of the achieved stratigraphic applications, several novel subjects can be discussed. First, the Python application seems to be an ergonomic and easy-access tool for visualizing the 3D stratigraphic architecture. The 3D models can be viewed in a web browser, so they are easily accessible to every user. The interactive 3D HTML models of essential stratigraphic elements (volumes and surfaces) allow for making quantitative measures, which is quite valuable to support many applied hydrogeological, geotechnical, and mining research interests. The second is the use of the granulometry as basic interpolating data (Figure 3), which is of special interest for geological explorations of porous sedimentary media because this physical parameter determines the storage of a fluid and its ability to move through the pores, including groundwater, oil, gas, brines, leachates, and pollutants. The third is related to a research gap concerning the need to program additional routines to assess the 3D models’ uncertainty, which is a task currently ongoing.

Regarding the Quaternary stratigraphic architecture and basement surface of the onshore LRD, additional findings to the above-described geological model have been obtained. For instance, the gravel and coarse sand lithosomes seem to be arranged into three stratigraphic intervals called lower (>50 m b.s.l.), middle (20–50 m b.s.l.), and upper (<20 m b.s.l.), as reflected in Figure 9. These coarse lithosomes seem to be deposited showing the following trend: in the distal LRD sector in the lower interval, in the central LRD sector in the middle interval, and in the entire LRD but with a great gravel lithosome in the embedding of the Llobregat River between the Garraf and Collserola massifs in the upper interval. This evolution fits well with the recent scientific literature (reflected in the lower cross-section in Figure 2) and seems to indicate a transgression during the lower and middle intervals with progressive retreating of coarse deposits towards proximal areas (getting closer and closer to the Garraf and Collserola massifs). This would be followed by a rapid regression with a spreading of coarse deposits toward distal areas to occupy the entire LRD. Figure 9 shows how the coarse deposits are not homogeneous in the entire LRD, as the old cylindric interpretation from the scientific literature indicated [11,35]. On the contrary, they fit very well with the recent interpretations, in which more parasequences with discontinuous coarse sedimentary bodies spread over the entire LRD were defined [7].

With respect to the basement, which is mostly composed of Pliocene rocks and Miocene and Paleozoic rocks to a minor extent [12,13,17], the different 3D views show a very irregular top surface. Figure 10 shows an overall steeped ‘horsts-grabens’ structure deepening toward the marine platform. This overall trend is interrupted with a raised (horst) sector in the El Prat de Llobregat town area, probably due to faulting, which is clearly visible in Figure 10C. In fact, when Figure 10C is compared to Figure 1, this horst sector narrowly coincides with the prolongation of the Montjuïc Mount elevation due to the action of the Morrot faults family below the LRD. The prolongation of the Tibidabo, Llobregat, and Morrot fault families could also explain the ‘horsts-grabens’ structure in the basement top surface of the onshore LRD in good agreement with recent research [7,47].

6. Conclusions

In 2004, the ACA created the METALL to compile and homogenize geological and hydrogeological information in order to assess the cumulative impact of the PDL civil infrastructures (with variable underground development) on the LRD groundwater bodies. On the basis of the high geological knowledge gained about the Quaternary onshore LRD, this area was selected to show the performance of a Python application that implements the packages NumPy, Pandas, Plotly, and Scipy for visualizing the 3D stratigraphic architecture. Although there are many commercial and open-source libraries and software devoted to geological visualization, its use often requires effort to learn how they work and money in the case of the commercial ones. The use of the Python language has proven efficient, ergonomic, and easy access to the results.

The Python application uses the granulometry dataset deduced from 433 borehole lithological records as interpolating data to create interactive 3D HTML models of essential stratigraphic elements. These models can be viewed in a web browser without any additional tools, so they are easily accessible to every user. The users can make different views, zoom, and move around, as well as hide any element by clicking that element in the right-side legend to see other details. The current development stage of the application enables three kinds of 3D HTML models presented in successive steps.

The first step models the 3D location of the boreholes granulometry (Figure 8). To this end, the georeferenced granulometry data were clustered into the three granulometry classes associated with three main lithosomes: gravel, coarse sand, and clay–silt. For better display, the main lithologies are colored blue for gravel, yellow for coarse sand, gray for silt–clay, and red for the basement, thus allowing for visual stratigraphic divisions by using conventional geological correlation procedures. The interactive 3D HTML model is included as supplementary material (3D_Boreholes_LRD.html).

The second step models the 3D volume of the gravel and coarse sand lithosomes arranged into three stratigraphic intervals, called lower (>50 m b.s.l.) in the distal LRD sector, middle (20–50 m b.s.l.) in the central LRD, and upper (<20 m b.s.l.) spread in the entire LRD (Figure 9). The main lithosome drawing uses the same colors for lithosomes. This 3D modeling fits well with the findings reported in the recent scientific literature. In general, a transgression during the lower and middle intervals followed by a rapid regression with a spreading of coarse sediments toward distal areas occupying the entire LRD is deduced. The interactive 3D HTML model is included as supplementary material (3D_Lithosomes_LRD.html).

The third step models the 3D basement (mostly Pliocene rocks and other older rocks to a minor extent) top surface of the Quaternary LRD formations (Figure 10). This model is created through the 87 boreholes reaching the basement. The mapped surface shows an overall steeped ‘horsts-grabens’ structure, deepening toward the marine platform and interrupted with a horst sector in the El Prat de Llobregat town area, probably due to faulting. This 3D modeling fits also well with the findings reported in the recent scientific literature. The interactive 3D HTML model is included as supplementary material (3D_Basement_LRD.html).

All that has been exposed about the ability to predict the volumes of the onshore Quaternary lithosomes and the basement top surface demonstrate the good performance of the Python programming language, and how the packages are a friendly software platform for visualizing the 3D stratigraphic architecture. The use of parameter granulometry as basic interpolating data widens the interest of this Python application in geological explorations of porous sedimentary media. Visualizing 3D models of essential stratigraphic elements is crucial in many applied hydrogeological, geotechnical, and mining research interests. The LRD is a porous sedimentary media forming a strategic, highly stressed groundwater body close to Barcelona city, so new geological applications of hydrogeological interest are especially welcome. For instance, the Python application may assist with crucial concerns for groundwater management and governance such as drilling of proper pumping wells, optimization of groundwater monitoring networks, as well as the evaluation of how the large civil works may impact the groundwater resource. The spatial distribution of granulometry may be of assistance to design compensatory measures for aquifer protection and recoveries such as artificial aquifer recharging and positive hydraulic barriers to reduce the advance of seawater intrusion and the mobilization of pollutants.

Currently, additional research on program routines aimed to assess the mapping uncertainty is ongoing. This paper contributes to the field of quantitative geological tools, which represent a poorly documented piece of information in the applied geology scientific literature.