GIS-Based Subsurface Analysis and 3D Geological Modeling as a Tool for Combined Conventional Mining and In-Situ Coal Conversion: The Case of Kardia Lignite Mine, Western Greece

: The development of three-dimensional geological models has proven to be critical for conceptualizing complex subsurface environments. This is crucial for mining areas due to their various hazards and unstable conditions. Furthermore, three-dimensional (3D) models can be the initial step for the development of numerical models in order to support critical decisions and sustainable mining planning. This paper illustrates the results and the development phases of a 3D geological model within the boundaries of the Kardia lignite deposit in western Macedonia, Greece. It also highlights the usefulness of a Geographic Information System (GIS) methodology in the subsurface geological and hydrogeological analysis regarding the Underground Coal Gasiﬁcation (UCG) methodology. In addition, the work focuses on the integrated geospatial framework that was developed to support the Coal-to-Liquids Supply Chain (CLSC) integration in unfavorable geological settings. A 3D subsurface geological model of the study area was developed to identify a suitable area for in situ coal conversion and UCG considering criteria related to speciﬁc coal thickness and depth. In this context, the suggested integrated geomodelling workﬂow can positively contribute to the implementation of conventional and innovative mining, saving time and reducing the cost to improve the quality of information needed to support decisions related to UCG implementation.


Introduction
Energy costs and demand have risen significantly during the recent COVID-19 and Ukraine crises (Figure 1) [1].According to the International Energy Agency (IEA) [2], in 2021, high-level increases in natural gas prices have prompted a significant switch to coal rather than natural gas to produce electricity in the energy markets of Europe, Asia, and the United States of America.Underground Coal Gasification (UCG) is a promising technology that converts coal into product gas, being more environmentally friendly than conventional mining techniques in subsurface environments.
Correctly managed UCG is an advanced clean coal technology that offers a potential solution for these challenges and broad application prospects.UCG opportunity arises in the energy transition period, where there is a need to reduce emissions significantly to comply with climate change commitments [3].UCG is not a new technology; it involves the gasification of unmined coal and the production of synthetic gas (syngas) for use in power generation or as chemical feedstock [4,5].UCG is appealing to expand recoverable coal reserves.It provides access to coal deposits that would otherwise remain unused due to complex and challenging mining and geological conditions (i.e., too deep/steep, low grade, thin seams) using conventional mining technologies [6,7].Correctly managed UCG is an advanced clean coal technology that offers a potential solution for these challenges and broad application prospects.UCG opportunity arises in the energy transition period, where there is a need to reduce emissions significantly to comply with climate change commitments [3].UCG is not a new technology; it involves the gasification of unmined coal and the production of synthetic gas (syngas) for use in power generation or as chemical feedstock [4,5].UCG is appealing to expand recoverable coal reserves.It provides access to coal deposits that would otherwise remain unused due to complex and challenging mining and geological conditions (i.e., too deep/steep, low grade, thin seams) using conventional mining technologies [6,7].
The main challenge is the development of competitive technologies for the production of synthesis gas and the production of electricity, heat, and synthetic liquid fuels on their basis [8].By combining UCG with Gas to Liquids (GTL) technologies, high-quality synthetic liquid transportation fuels can produce environmental and economic advantages over alternate conventional and unconventional processes [9].
Despite the extensive research, there are many areas of improvement of the UCG process and adaptation to local conditions that must be addressed through additional research.These improvements need to advance the selection of the best target area and the effectiveness of the gasification process, making sure it does not impact the environment negatively.
Many countries have already researched implementing UCG technologies, including the long-term commercial operation of several UCG plants in the former Soviet Union.However, there is limited information on applying this technology in Greece.In the ODYSSEUS project (Coal-to-liquids supply chain integration because of operational, economic, and environmental risk assessments under unfavorable geological settings), we have identified and reviewed the criteria for selecting potential UCG locations.
Studying a subsurface environment within the boundaries of a large area is a complex task.This procedure consists of many parameters, which can be examined separately under conventional methods but needs added workload and resources.On the other hand, a developed relational database can manipulate and manage various data related to GIS approaches.Furthermore, recent publications have shown that integrating a geospatial framework is an efficient tool for analyzing and mapping subsurface geological and hydrogeological environments [10,11].
In addition, geomodelling requires GIS databases and tools [12], due to their ability to interact with the data spatially, to understand the subsurface and its mechanisms better, and visualize it in 2D and 3D dimensions [11,13].Even though much work has been done in the last century concerning geomodelling and developed GIS databases, there is no standard widely accepted workflow.Due to the high cost of various software packages The main challenge is the development of competitive technologies for the production of synthesis gas and the production of electricity, heat, and synthetic liquid fuels on their basis [8].By combining UCG with Gas to Liquids (GTL) technologies, high-quality synthetic liquid transportation fuels can produce environmental and economic advantages over alternate conventional and unconventional processes [9].
Despite the extensive research, there are many areas of improvement of the UCG process and adaptation to local conditions that must be addressed through additional research.These improvements need to advance the selection of the best target area and the effectiveness of the gasification process, making sure it does not impact the environment negatively.
Many countries have already researched implementing UCG technologies, including the long-term commercial operation of several UCG plants in the former Soviet Union.However, there is limited information on applying this technology in Greece.In the ODYSSEUS project (Coal-to-liquids supply chain integration because of operational, economic, and environmental risk assessments under unfavorable geological settings), we have identified and reviewed the criteria for selecting potential UCG locations.
Studying a subsurface environment within the boundaries of a large area is a complex task.This procedure consists of many parameters, which can be examined separately under conventional methods but needs added workload and resources.On the other hand, a developed relational database can manipulate and manage various data related to GIS approaches.Furthermore, recent publications have shown that integrating a geospatial framework is an efficient tool for analyzing and mapping subsurface geological and hydrogeological environments [10,11].
In addition, geomodelling requires GIS databases and tools [12], due to their ability to interact with the data spatially, to understand the subsurface and its mechanisms better, and visualize it in 2D and 3D dimensions [11,13].Even though much work has been done in the last century concerning geomodelling and developed GIS databases, there is no standard widely accepted workflow.Due to the high cost of various software packages and lack of data and resources, 3D geological modeling of the subsurface has not been applied frequently in publications and various studies [14,15].
Within such a framework, the present study aims to develop an efficient methodology of integrating borehole and various data to map and research the spatial distribution of lignite seam thickness in the Kardia mine.The developed 3D geological model was integrated with specific geospatial analysis tools for the lignite deposit and was applied to locate and identify potentially suitable sites for UCG implementation employing combined conventional mining and in situ coal conversion.

Geological and Geomorphological Settings
The lignite mine is located in the Ptolemais Lignite Basin, dominated by E-W trending normal Quaternary faults [16].The Ptolemais Basin covers a surface area of approximately 600 km 2 .The longitudinal axis of the lignite segment of the Ptolemais Basin, with an NW-SE direction, exceeds 20 km in length, while the maximum width reaches about 20 km.The basin (Figure 2) is filled with Late Miocene to Pleistocene lake sediments, including intercalated lignites and alluvial deposits with a total thickness of up to 600 m.
The sediments consist of sandy clays, sandy marls, clayey marls, calcareous sands, and conglomerates.Above lay Quaternary conglomerates of terrestrial and fluvioterrestrial origin.The exposed stratigraphic sequence in the Kardia lignite mine belongs to the Early Pliocene and Ptolemais Formation.Based on the borehole data, the stratigraphic sequence in the lignite mine has an average thickness of about 300 m and consists of an alternation of lignite and marls, with intercalated sands and silts [17][18][19].The Ptolemais Basin is part of a tectonic trench over 250 km extending from northern Greece into North Macedonia [20].
Within such a framework, the present study aims to develop an efficient methodology of integrating borehole and various data to map and research the spatial distribution of lignite seam thickness in the Kardia mine.The developed 3D geological model was integrated with specific geospatial analysis tools for the lignite deposit and was applied to locate and identify potentially suitable sites for UCG implementation employing combined conventional mining and in situ coal conversion.

Geological and Geomorphological Settings
The lignite mine is located in the Ptolemais Lignite Basin, dominated by E-W trending normal Quaternary faults [16].The Ptolemais Basin covers a surface area of approximately 600 km 2 .The longitudinal axis of the lignite segment of the Ptolemais Basin, with an NW-SE direction, exceeds 20 km in length, while the maximum width reaches about 20 km.The basin (Figure 2) is filled with Late Miocene to Pleistocene lake sediments, including intercalated lignites and alluvial deposits with a total thickness of up to 600 m.
The sediments consist of sandy clays, sandy marls, clayey marls, calcareous sands, and conglomerates.Above lay Quaternary conglomerates of terrestrial and fluvioterrestrial origin.The exposed stratigraphic sequence in the Kardia lignite mine belongs to the Early Pliocene and Ptolemais Formation.Based on the borehole data, the stratigraphic sequence in the lignite mine has an average thickness of about 300 m and consists of an alternation of lignite and marls, with intercalated sands and silts [17][18][19].The Ptolemais Basin is part of a tectonic trench over 250 km extending from northern Greece into North Macedonia [20].According to previous literature, the thickness of the lignite-bearing layers (including intercalations) in the adjacent mining area ranges between 80 m and 140 m at the western boundary of the mine near the "Askio" (Siniatsiko) mountain.After that, it increases to the SW to more than 210 m, where the overlying strata have a thickness of about 150 m.In the central and north-western parts of the mine, the thickness of the overlying strata is about 20 m-60 m, and the thickness of the lignite seams varies [22].
Regarding the Digital Elevation Model (DEM), the Kardia mine is an open and excavated area with surface elevations ranging between +622 m and +812 m above sea level (a.s.l.) from E to W direction, respectively.DEM is the most common digital 3D representation of a topographic surface concerning any reference datum.Therefore, DEM is frequently used to determine terrain attributes that include elevation at any point, slope, Mining 2022, 2 300 and aspect [23].In our study area, DEM was imported from satellite data using the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), with a nominal horizontal accuracy of 15 m and 20 m, respectively [24].
In addition, the surface elevation value for each borehole was obtained from the DEM dataset.The mining area covers approximately 18 km 2 near the Pontokomi village (Figure 3).The lignite basin of the deposit progressively broadens from W to E with a gradational decrease in surface altitudes.
In the central and north-western parts of the mine, the thickness of the overlying strata is about 20 m-60 m, and the thickness of the lignite seams varies [22].
Regarding the Digital Elevation Model (DEM), the Kardia mine is an open and excavated area with surface elevations ranging between +622 m and +812 m above sea level (a.s.l.) from E to W direction, respectively.DEM is the most common digital 3D representation of a topographic surface concerning any reference datum.Therefore, DEM is frequently used to determine terrain attributes that include elevation at any point, slope, and aspect [23].In our study area, DEM was imported from satellite data using the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), with a nominal horizontal accuracy of 15 m and 20 m, respectively [24].
In addition, the surface elevation value for each borehole was obtained from the DEM dataset.The mining area covers approximately 18 km 2 near the Pontokomi village (Figure 3).The lignite basin of the deposit progressively broadens from W to E with a gradational decrease in surface altitudes.

Geospatial Analysis Framework and Resolution
Rockware's Rockworks17 and ESRI ArcGIS 10.4 software were used for geomodelling and spatial analysis processing tasks.The Rockworks software is a geomodeling software able to create block models and cross-sections in an area, providing the option to convert the exported data into ArcGIS shape easily or grid files and vice versa as input data.Thus, for the 3D geological modeling in the Rockworks environment, the lateral

Geospatial Analysis Framework and Resolution
Rockware's Rockworks17 and ESRI ArcGIS 10.4 software were used for geomodelling and spatial analysis processing tasks.The Rockworks software is a geomodeling software able to create block models and cross-sections in an area, providing the option to convert the exported data into ArcGIS shape easily or grid files and vice versa as input data.Thus, for the 3D geological modeling in the Rockworks environment, the lateral blending technique or Lithoblending algorithm was implemented for interpolating geological units' spatial distribution regarding horizontal and vertical (values) axes [25].
In this study, two spatial data categories have been utilized: (i) discrete representation of objects and (ii) representation with continuous fields.The type of discrete representation consists of well-defined objects in terms of shape, color, and form.Under this context, a geographical area can be described as a region with lakes, rivers, faults, land cover, or geological formations representing polylines, points, or polygons.
In order to produce valuable and comparable outcomes, preparation and filtering of the original geospatial data were applied before the thorough spatial analysis.According to the following schematic diagram (Figure 4), a specific workflow was implemented hierarchically using spatial analysis and specific geoprocessing tools within the environment of ArcGIS 10.4 version (left part of Figure 4) [26] and Rockworks (right part of Figure 4).geological formations representing polylines, points, or polygons.
In order to produce valuable and comparable outcomes, preparation and filtering of the original geospatial data were applied before the thorough spatial analysis.According to the following schematic diagram (Figure 4), a specific workflow was implemented hierarchically using spatial analysis and specific geoprocessing tools within the environment of ArcGIS 10.4 version (left part of Figure 4) [26] and Rockworks (right part of Figure 4).According to the upper part of the workflow, the first step was the selection and correction of the Digital Elevation Model of the study area.Then, the DEM correction was applied in the ArcGIS environment using the Fill geoprocessing tool.According to ESRI, peaks and sinks are often errors due to the resolution of the data or rounding of elevations to the nearest integer value.The fill tool uses the equivalents of several tools, such as Focal Flow, Flow Direction, Sink, Watershed, and Zonal Fill, to locate and fill the sinks [27][28][29][30].Thus, the tool iterates until all peaks or sinks within the specified z limit are filled.Subsequently, after implementing the sink tool, more than 15 sink areas were identified with a total area of more than 15,000 m 2 and removed from the initial Digital Elevation Model dataset.After the correction with the Fill algorithm of the DEM, the elevation values were extracted to add the information of each borehole within the boundaries of the study area.According to the presented workflow, georeferencing and digitization of faults according to previous studies and available references were applied to be analyzed in Rockworks software as linear lineaments to produce rose diagrams (lower part of the flow chart).After the integration of borehole data in the geodatabase, followed the visualization and spatial analysis of the 304 boreholes dataset under the process of the study area in order to compare the observed point patterns to the ones generated by an independent random process (IRP), also known as complete spatial randomness (CSR).This technique was implemented using the Average Nearest Neighbour (ANN) tool from ArcGIS 10.4 version.According to the reference ANN algorithm, which is included in Arcmap software [31,32], checks whether or not the observed first-order nearest neighbor is consistent with a distribution of points if this process was completely random.The average nearest neighbor (ANN) tool computes the first nearest neighbor mean distance for all points [33,34].The following equation ( 1) is used to calculate the ANN ratio : where OMD is the observed mean distance between each feature and its nearest neighbor (2): and EMD (3) is the Expected Mean Distance for the features given in random pattern: In the above equations, d i equals the distance between feature i and its nearest neighboring feature, n corresponds to the total number of features, and A is the area of the minimum enclosing rectangle around all features or it can be specified.
The ANN z-score is calculated as (4): where SE is calculated as follows ( 5): ANN results showed that the Observed Mean Distance was approximately 134.4 m, and the Expected Mean Distance (EMD) was calculated at approximately 100 m.Also, the ANN ratio or Nearest Neighbor Ratio (NNR) was calculated as 1,355, and regarding the z-score (Figure 5), the spatial distribution pattern was categorized as dispersed.EMD is a key aspect of calculating the distance at which locations between boreholes (points) are no longer related where the spatial dependence ceases.The EMD as 100 m was considered to evaluate the horizontal and vertical distance spacing in Rockworks software and, therefore, develop the 3D geological model to have an efficient, accurate, but minimum computational effort discretization.According to the upper part of the workflow, a critical step was developing and homogenizing a catalog The EMD as 100 m was considered to evaluate the horizontal and vertical distance spacing in Rockworks software and, therefore, develop the 3D geological model to have an efficient, accurate, but minimum computational effort discretization.According to the upper part of the workflow, a critical step was developing and homogenizing a catalog with 304 boreholes, containing their exact location, elevation above sea level, total depth, and the lithological unit at each depth.
The description of the lithological units that were used and their codes are presented in Table 1.Nine lithological units were coded and visualized (Table 1).In addition, the Rockworks Borehole Manager tool was utilized to create the grid models for the subsurface of each lithological unit in order to visualize them in 2D cross-sections and 3D diagrams (lowest part of Figure 4).For the creation of the block model, in terms of the upper surface and lower subsurface filters, Inverse Distance Gridding (IDG) was implemented, which is suggested in Rockworks software and is based on Newton's Gravity Equation (6).According to this method, each value assigned to a grid node is a weighted average of all data points or several directionally distributed neighboring control points.Specifically, the value of each data point is weighted according to its inverse distance from the grid node, taken to a user-selected power.The greater the exponent specified value, the more localized the gridding since distant points will have less influence on the value assigned to each grid node.
Z-node is the value of the assigned node values, where the value of each data point is weighted according to the inverse of its distance d from the grid node, taken to the selected power (6).When the distance factor d is greater than one, the values of distant points have less influence than nearby points.For the upper and lower surface of the developed block model, the automatic methodology was selected based on borehole collar elevations, where the software automatically interpolated a surface representing the elevations at the top of the boreholes (using current project dimensions and the inverse-distance gridding method).
The advantage of using the IDG approach is that it can produce a smooth and continuous grid without exaggeration regarding extrapolations beyond the given data points.On the other hand, IDG can create a "bulls-eye" effect in some areas of the produced dataset.According to Rockware, the range of grid values will be smaller than the data point range.In other words, the highest value of the grid will be less than the maximum data point, and the lowest grid value will be greater than the minimum data point [25].
Regarding the workflow (Figure 4), statistical analysis has been applied to filter the suitable boreholes according to the UCG criteria.The primary criteria that have been used refer to the coal seam thickness, where single seam thickness should be at least 1 m [35,36] and coal seam depth greater than 200 m to minimize potential gases [37,38].After filtering the 304 boreholes, 10 boreholes in a specific region were selected, where a new high-resolution lithological model was produced to illustrate coal seams better and visualize the suggested subsurface environment in greater detail.
The implemented methodology provides a distinct and integrated workflow regarding how the lithological model has been developed in a specific area.The approach followed regarding Greek open-pit lignite regions needs specific steps to be implemented to analyze the geological subsurface in the context of the UCG technique.In addition, this workflow attempted to present a specific flow in terms of how geomodelling combines GIS databases and tools to understand better the subsurface and its mechanisms, illustrating it in 2D and 3D dimensions according to UCG criteria.Previous works have been published in lithological modeling concerning geomodelling and developed GIS databases, but there is no standard widely accepted workflow [14].
Thus, due to the high cost of various software packages and lack of data and resources, 3D geological modeling of the subsurface has not been widely utilized in publications and various studies [15].

Geospatial Analysis Results
An essential part of the work was the investigation and the georeferencing process of the faults' traces (red polylines) on the surface of the mine area, based on the literature and geological studies.Most of the tectonic lines within the boundaries of the studied area were digitized and imported into the database.More specifically, more than 300 faults were mapped (Figure 6) in the broader area, where 31 of them were inside the polygon of the mine.According to the georeferenced faults' traces (Figure 6) and the following rose diagram (Figure 7), faults dissect the mining field into several fault blocks.According to the georeferenced faults' traces (Figure 6) and the following rose diagram (Figure 7), faults dissect the mining field into several fault blocks.According to the georeferenced faults' traces (Figure 6) and the following r gram (Figure 7), faults dissect the mining field into several fault blocks.The highest percentage of all faults' lineation (>150 faults segments) by length is following an NW-SE direction in the broader area of the mine.The total length of the faults is 224 km, the average bin population is 51.81, and the standard deviation of the bin population is 30.22.

Three-Dimensional Geological Modelling
Currently, 3D geological modeling is a rapidly expanding discipline as demand for knowledge about the structure and composition of the subsurface is increasing.In order to create the 3D model, existing borehole data were grouped in lithological units to be imported into the RockWorks Borehole Manager.
After grouping the lithological layers, the study area boundaries were calculated with the Minimum/Maximum Eastern and Northern limits.Then, based on the borehole data and the boreholes' coordinates, a grid model was created with specific dimensions and nodes.The Spatial Reference System utilized was the Universal Transverse Mercator (UTM) 34 (Northern Hemisphere)/WGS 84.It is worth mentioning that lithological units are repeated with depth in each borehole (sand, clay, coal, or other formations).Therefore, the following procedure coded the lithological types like integer or real number values, based on the "G" value for each rock/soil type.For instance, "CO = coal" was coded with a "2" and "sand" with a "4.5".The "G" value is just a numeric value that Rockworks uses to represent the material in lithology models.A sample 3D profile of the T6D-070 borehole illustrates the lithological layers using a RockPlot3-D window in the following image (Figure 8).
The three-dimensional lithological model (Figure 9) was generated taking into account three primary characteristics of the solids: extent, base, and top-level of each lithological layer, depth, thickness, and descriptive data related to lithological type.For the better interpolation of the model, a gradational method was applied to standardize and group similar lithological types.The 3D block model represents the different lithology types, creating a 3D dimensional voxel diagram that illustrates these lithological types.Lithology models are displayed as all-voxel diagrams, with each color-coded voxel based on its assigned material type.For example, in our project, the black color corresponds to a coal formation and the yellow color to a marl formation.
are repeated with depth in each borehole (sand, clay, coal, or other formations).Therefore, the following procedure coded the lithological types like integer or real number values, based on the "G" value for each rock/soil type.For instance, "CO = coal" was coded with a "2" and "sand" with a "4.5".The "G" value is just a numeric value that Rockworks uses to represent the material in lithology models.A sample 3D profile of the T6D-070 borehole illustrates the lithological layers using a RockPlot3-D window in the following image (Figure 8).The three-dimensional lithological model (Figure 9) was generated taking into account three primary characteristics of the solids: extent, base, and top-level of each lithological layer, depth, thickness, and descriptive data related to lithological type.For the better interpolation of the model, a gradational method was applied to standardize and group similar lithological types.The 3D block model represents the different lithology types, creating a 3D dimensional voxel diagram that illustrates these lithological types.Lithology models are displayed as all-voxel diagrams, with each color-coded voxel based on its assigned material type.For example, in our project, the black color corresponds to a coal formation and the yellow color to a marl formation.
The development of the model was calculated in three dimensions as a 3D voxel/solid diagram.A 3D model of lithological variations within the basin was developed by extrapolating 304 boreholes using a 3-dimensional gridding process.Subsurface lithology was defined through the identification of distinctive lithological packages.A log correlation was performed using the lithoblending algorithm to create the model.This unique algorithm exports lithology types outward from the borehole data and trenches into a solid The development of the model was calculated in three dimensions as a 3D voxel/solid diagram.A 3D model of lithological variations within the basin was developed by extrapolating 304 boreholes using a 3-dimensional gridding process.Subsurface lithology was defined through the identification of distinctive lithological packages.A log correlation was performed using the lithoblending algorithm to create the model.This unique algorithm exports lithology types outward from the borehole data and trenches into a solid block model.This approach produces "lithozones" around each borehole which can end suddenly when a lithological zone from a near borehole is encountered [39].More specifically, this algorithm acts as a nearest neighbor geospatial methodology, exclusively designed for lithological modeling in the Rockworks environment.The lithoblending algorithm detects the horizontal plane around each log.Then, it iteratively assigns the corresponding numeric value, which has been set in the lithology table, to the next node unless it bumps into a node that has already been interpolated, following a growing circle around each log.It should be noted that the horizontal weighting capacity of the lithoblending algorithm activated the illustration of the horizontal lenses from the different formations as frequently observed in environments or basins with fluvial sedimentation [40].
Available subsurface data provided sufficient detail to distinguish major lithological types confidently and with sufficient detail within these units to develop a subsurface geologic model.In our case (Figure 9), 9 primary lithological layers were depicted across the study area.In particular, breccia, conglomerate, sand/sandstone, clay, marl, coal, limestone, soil, and silt/siltstone formations were displayed in the 3D block model.In the study area, 10 vertical and horizontal cross-sections were produced (Figure 9) for the 3D model displaying the results of the numerical interpolation as a fence diagram that highlights lithological variations within the basin fill.
The cross sections' colors correspond to the lithological units in the fence diagram (Figure 10).The upper surfaces of the cross-sections were then visualized according to the elevation of each borehole (highest 732.96 m and lowest 634.13 m).Finally, the model base was visualized according to the borehole data's lowest elevations (the difference between the depth and elevation).According to the derived model, the marly, sandy, clayey formations and coal layers seem to dominate the basin's subsurface, according to the cross-sections.
rithm detects the horizontal plane around each log.Then, it iteratively assigns the corresponding numeric value, which has been set in the lithology table, to the next node unless it bumps into a node that has already been interpolated, following a growing circle around each log.It should be noted that the horizontal weighting capacity of the lithoblending algorithm activated the illustration of the horizontal lenses from the different formations as frequently observed in environments or basins with fluvial sedimentation [40].Available subsurface data provided sufficient detail to distinguish major lithological types confidently and with sufficient detail within these units to develop a subsurface geologic model.In our case (Figure 9), 9 primary lithological layers were depicted across the study area.In particular, breccia, conglomerate, sand/sandstone, clay, marl, coal, limestone, soil, and silt/siltstone formations were displayed in the 3D block model.In the study area, 10 vertical and horizontal cross-sections were produced (Figure 9) for the 3D model displaying the results of the numerical interpolation as a fence diagram that highlights lithological variations within the basin fill.
The cross sections' colors correspond to the lithological units in the fence diagram (Figure 10).The upper surfaces of the cross-sections were then visualized according to the elevation of each borehole (highest 732.96 m and lowest 634.13 m).Finally, the model base was visualized according to the borehole data's lowest elevations (the difference between the depth and elevation).According to the derived model, the marly, sandy, clayey formations and coal layers seem to dominate the basin's subsurface, according to the crosssections.According to the previous sections, specific tools and methodologies have been used to create the 3D geological model and manage the geospatial data in the study area.The result of this integration was the geological interpretation and modeling of the study area's subsurface environment within the boundaries of the Kardia lignite mine regarding the 304 borehole data and its spatial distribution.The relational database, which has been According to the previous sections, specific tools and methodologies have been used to create the 3D geological model and manage the geospatial data in the study area.The result of this integration was the geological interpretation and modeling of the study area's subsurface environment within the boundaries of the Kardia lignite mine regarding the 304 borehole data and its spatial distribution.The relational database, which has been developed, contains the digitized faults traces, borehole data, 2D and 3D maps, and datasets for future analysis to enhance the geological and hydrogeological understanding of the area under examination.

UCG Site Identification
The next step was to identify and delineate a proposed subsurface suitable area for combined conventional mining and in situ coal conversion within the boundaries of the created block lithological model, consisting of 304 borehole data.According to reference and UCG technology, two significant criteria were considered (i) borehole location with a total depth greater than 200 m and (ii) coal thickness greater than 1 m.A series of Structured Query Language (SQL) queries combined with statistical results were implemented into Rockworks in order to filter the 304 boreholes to find suitable areas fitting the UCG abovementioned criteria, considering the needed depth and coal continuity.After the integrated analysis and the spatial distribution of the different lithologies in compliance with the primary criteria, a particular region was identified consisting of only 10 boreholes (Table 2).The developed lithological model was filtered and cut off in order to surround these borehole locations and their lithological distribution.Specifically, the developed block model was isolated with spatial filtering to contain the 10 boreholes and reconstructed in higher resolution, considering horizontal and vertical distances for better visualization and interpretation.In this framework, the data of the 10 boreholes (Table 2): KND-5/14, KND-3/14, T6D-040, T6D-134, T6D-050, 242/228, T6-D60, T6D-059, T6D-051, and T6D-079 (Table 3), were considered, incorporating the interpolated lithological sequence in the suggested area, and the alternation of lignite and marls, with intercalated sands and silts.The results validated and highlighted more clearly the sequence corresponding to the Late Miocene to Pleistocene formations, which fill the basin, including intercalated lignite and recent deposits with a total thickness of up to 600 m.Two distinct lignite horizons have been identified (Figure 11), altering with marly and clayey formations.The developed new high-resolution solid model and the obtained cross-section from W to E (Figure 12) were generated with 1 m vertical and horizontal distances as long as 10 m horizontal spacing.
According to the following cross-section from W to E direction (Figure 12), a specific area with lignite lithological units (black dashed line) was delineated and selected with an extent of approximately 135,000 m 2 .Top and bottom elevation heights for the volumetric of the coal were set from 500 m to 350 m a.s.l.(185 m to 325 m depth), respectively.In addition, spatial filtering was applied in X-Y coordinates from 564,600 m to 565,500 m and from 4,474,450 m to 4,474,600 m (UTM 34 N).
and the alternation of lignite and marls, with intercalated sands and silts.The results validated and highlighted more clearly the sequence corresponding to the Late Miocene to Pleistocene formations, which fill the basin, including intercalated lignite and recent deposits with a total thickness of up to 600 m.Two distinct lignite horizons have been identified (Figure 11), altering with marly and clayey formations.The developed new highresolution solid model and the obtained cross-section from W to E (Figure 12) were generated with 1 m vertical and horizontal distances as long as 10 m horizontal spacing.According to the following cross-section from W to E direction (Figure 12), a specific area with lignite lithological units (black dashed line) was delineated and selected with an extent of approximately 135,000 m 2 .Top and bottom elevation heights for the volumetric of the coal were set from 500 m to 350 m a.s.l.(185 m to 325 m depth), respectively.In addition, spatial filtering was applied in X-Y coordinates from 564,600 m to 565,500 m and from 4,474,450 m to 4,474,600 m (UTM 34 N).Based on the created cross-section of the solid model (Figure 12), the total volume of the coal in the suggested area (black dashed outline) was calculated as more than 4 * 10 6 m 3 .In addition, Figure 13 illustrates a 3D diagram in the same area (red rectangle), visualizing only the isolated and selected lignite layers of the developed high-resolution block model with vertical exaggeration x6.Based on the created cross-section of the solid model (Figure 12), the total volume of the coal in the suggested area (black dashed outline) was calculated as more than 4 * 10 6 m 3 .In addition, Figure 13   Based on the created cross-section of the solid model (Figure 12), the total volume of the coal in the suggested area (black dashed outline) was calculated as more than 4 * 10 6 m 3 .In addition, Figure 13 illustrates a 3D diagram in the same area (red rectangle), visualizing only the isolated and selected lignite layers of the developed high-resolution block model with vertical exaggeration x6.

Three-Dimensional Hydrolithological Model
The main aquifer system is developed in the sand, conglomerate, and clay alternations consisting of numerous confining layers that are delimited by impermeable clay with the piezometric level ranging from 640 m to 660 m a.s.l.[42].

Three-Dimensional Hydrolithological Model
The main aquifer system is developed in the sand, conglomerate, and clay alternations consisting of numerous confining layers that are delimited by impermeable clay with the piezometric level ranging from 640 m to 660 m a.s.l.[42].
According to pumping test results, in the broader area of Ptolemais mines, the hydraulic conductivity of clay and sand alternations ranges from 1.03 to 25.92 m/day, and marl and lignite alternations range from 0.05 to 0.32 m/day (Table 3) [42,43].The main hydrolithological units, which are the Pleistocene-Pliocene formation and the lignite sequence, were identified based on the 3D lithological block models in the targeted site.Within such a framework, the 3D model of the main hydrolithological units was visualized based on the 3D lithological model by matching each hydrolithological unit with its respective average hydraulic conductivity (Figure 14).
According to the following hydrolithological cross-section (Figure 15), the maximum depth of the main aquifer developed in Pleistocene-Pliocene formations is approximately 150 m from the surface.geted site.Within such a framework, the 3D model of the main hydrolithological unit was visualized based on the 3D lithological model by matching each hydrolithologica unit with its respective average hydraulic conductivity (Figure 14).
According to the following hydrolithological cross-section (Figure 15), the maximum depth of the main aquifer developed in Pleistocene-Pliocene formations is approximately 150 m from the surface.

Conclusions
It is well known that climate policies are becoming increasingly stringent, and closure phase of lignite mines has already started in Europe.In addition, regions adja to lignite mines face an unpredictable future.Therefore, the transition phase mus based on proper planning that considers the existing infrastructure and social awaren The main objective of this study is to highlight the geospatial capabilities for geolog and hydrogeological conceptualization within the boundaries of the Greek lignite min Kardia regarding the UCG approach.
The development of a geological model plays a key role in better understanding subsurface regarding the study area.In addition, the created database can offer hel information for future investigation and technical studies utilizing the existing spatia tabase and the produced datasets.Thus, 3D geological models illustrate the import of 3D subsurface modeling in a complex geological environment underlain by multi-l lignite horizons.Spatial results of the developed lithological model revealed the most able location for combined conventional mining and in situ coal conversion in the so

Conclusions
It is well known that climate policies are becoming increasingly stringent, and the closure phase of lignite mines has already started in Europe.In addition, regions adjacent to lignite mines face an unpredictable future.Therefore, the transition phase must be based on proper planning that considers the existing infrastructure and social awareness.The main objective of this study is to highlight the geospatial capabilities for geological and hydrogeological conceptualization within the boundaries of the Greek lignite mine of Kardia regarding the UCG approach.
The development of a geological model plays a key role in better understanding the subsurface regarding the study area.In addition, the created database can offer helpful information for future investigation and technical studies utilizing the existing spatial database and the produced datasets.Thus, 3D geological models illustrate the importance of 3D subsurface modeling in a complex geological environment underlain by multi-layer lignite horizons.Spatial results of the developed lithological model revealed the most suitable location for combined conventional mining and in situ coal conversion in the southwest region of the mine.
This work also proves that developing a 3D modeling approach is an efficient way to capture the subsurface complexity of unfavorable geological settings in the context of an integrated approach and better knowledge of the subsurface regime.Hence, identifying the most suitable region for UCG methodology according to the two major criteria, as long as creating a relational database with useful datasets, can be utilized in future studies and other research projects.A limitation of this work is that UCG methodology is more complex regarding its criteria and needs further study and geotechnical analysis within the boundaries of the suggested region.Although the developed high-resolution model, which filtered down the initial block diagram, seems to be an efficient approach to map the various lithological units and their spatial distribution.In this way, they can support sustainable practices and policies within and around the boundaries of lignite mines.To summarize, the methodology that was applied in the study area consists of a fivefold process:

•
The creation of a geoinventory and the development of a GIS-managed database;

•
The generation of 2D maps and 3D models as long as datasets;

•
The creation of 3D lithological and hydrolithological models through a solid modeling methodology and regarding spatial analysis and SQL queries; The identification of a suitable area regarding borehole information, UCG major criteria, and spatial distribution of lithology,

•
The development of a new high-resolution model to better interpret and understand the subsurface conditions.
Further geological and geotechnical studies can be implemented using the existing relational database and the developed 3D lithological and hydrolithological models regarding UCG methodology.Furthermore, this work showed that the iterated combination of GIS and Rockworks algorithms could be used efficiently to understand better and interpret complex subsurface parameters and mechanisms, providing a key role for UCG technology and local decision-making.From this point of view, 3D modeling seems to be the most satisfying way to interpret the subsurface complexity of the geological environment.In addition, the generated block models can be used in future studies of thorough groundwater flow and contaminant transport, where the developed databases could be used as required input data.
In conclusion, modern approaches and technologies can be combined to be transformed into diversified and sustainable resources that provide a better and more sustainable future regarding the environmental, economic, and public aspects.Therefore, the presented workflow can be applied to develop conceptual and initial models in sedimentary basins and complex geological environments, illustrating a meaningful and helpful approach combining results from conceptual models and GIS datasets.

Figure 1 .
Figure 1.Crude oil prices from 2013 to 2022, based on data from US Energy Information Administration-EIA [1].COVID-19 and recent Ukraine crisis-derived price increases are shown (red arrows).

Figure 1 .
Figure 1.Crude oil prices from 2013 to 2022, based on data from US Energy Information Administration-EIA [1].COVID-19 and recent Ukraine crisis-derived price increases are shown (red arrows).

Figure 2 .
Figure 2. Simplified geological map of the Ptolemais area, located NW and SE of the city of Ptolemais (based and modified on Institute of Geology and Mineral Exploration (IGME) geological sheet, Ptolemais 1:50.000,[21]).

Figure 2 .
Figure 2. Simplified geological map of the Ptolemais area, located NW and SE of the city of Ptolemais (based and modified on Institute of Geology and Mineral Exploration (IGME) geological sheet, Ptolemais 1:50.000,[21]).

Figure 3 .
Figure 3.The map represents the spatial location of the analyzed boreholes (304) combined with the Digital Elevation Model in the study area using ArcGIS 10.4 version.The base map is provided by Sentinel 2b satellite European Union-Copernicus program) with sensing date: 31 August 2021.

Figure 3 .
Figure 3.The map represents the spatial location of the analyzed boreholes (304) combined with the Digital Elevation Model in the study area using ArcGIS 10.4 version.The base map is provided by Sentinel 2b satellite European Union-Copernicus program) with sensing date: 31 August 2021.

Figure 4 .
Figure 4. Workflow of the data preparation for the geospatial analysis and creation of the database using various datasets.According to the upper part of the workflow, the first step was the selection and correction of the Digital Elevation Model of the study area.Then, the DEM correction was applied in the ArcGIS environment using the Fill geoprocessing tool.According to ESRI, peaks and sinks are often errors due to the resolution of the data or rounding of elevations to the nearest integer value.The fill tool uses the equivalents of several tools, such as Focal Flow, Flow Direction, Sink, Watershed, and Zonal Fill, to locate and fill the sinks[27][28][29][30].

Figure 4 .
Figure 4. Workflow of the data preparation for the geospatial analysis and creation of the database using various datasets.

Figure 5 .
Figure 5. Distribution pattern of boreholes regarding Average Nearest Neighbor (ANN) report.

Figure 5 .
Figure 5. Distribution pattern of boreholes regarding Average Nearest Neighbor (ANN) report.

Mining 2022, 2 ,Figure 6 .
Figure 6.The map represents the spatial location of the digitized faults in the broader area of the Kardia lignite mine.The base map is provided by Sentinel 2b satellite European Union-Copernicus program) with sensing date: 31 August 2021.

Figure 6 .
Figure 6.The map represents the spatial location of the digitized faults in the broader area of the Kardia lignite mine.The base map is provided by Sentinel 2b satellite European Union-Copernicus program) with sensing date: 31 August 2021.

Figure 6 .
Figure 6.The map represents the spatial location of the digitized faults in the broader ar Kardia lignite mine.The base map is provided by Sentinel 2b satellite European Union-Co program) with sensing date: 31 August 2021.

Figure 7 .
Figure 7. Rose diagram showing the distribution of fault trends at the broader area of Kar mine.

Figure 7 .
Figure 7. Rose diagram showing the distribution of fault trends at the broader area of Kardia mine.

Figure 8 .
Figure 8. Sample profile of the T6D-070 borehole (Vertical Exaggeration =3x, 1350/320 compass bearing/angle from the horizon).The development of the 3-D model of the study area was created with RockWorks version 17 (standard license).

Figure 8 .
Figure 8. Sample profile of the T6D-070 borehole (Vertical Exaggeration =3x, 1350/320 compass bearing/angle from the horizon).The development of the 3-D model of the study area was created with RockWorks version 17 (standard license).

Figure 9 .
Figure 9. Interpolation of the 3D block model (Vertical Exaggeration = 8x, 1100/260 compass bearing/angle from the horizon).The development of the 3-D model of the study area was created with RockWorks version 17 [41].Vertical axis represent a.s.l values.

Figure 10 .
Figure 10.Grid lithology fence panels (Horizontally and vertically; 262 • /400) are displayed in 3D using the background colors based on the lithology patterns.

Figure 11 .
Figure 11.High-resolution 3D lithological block model in the suggested area according to UCG criteria.Vertical axis represent a.s.l values.

Figure 11 . 14 Figure 12 .
Figure 11.High-resolution 3D lithological block model in the suggested area according to UCG criteria.Vertical axis represent a.s.l values.Mining 2022, 2, FOR PEER REVIEW 14

Figure 12 .
Figure 12.W-E cross-section on selected boreholes with Vertical Exaggeration x1.7.The black dashed line represents the delineated lignite area for the UCG approach.Vertical axis represent elevation above sea level (a.s.l) and depth values.
illustrates a 3D diagram in the same area (red rectangle), visualizing only the isolated and selected lignite layers of the developed high-resolution block model with vertical exaggeration x6.
elevation above sea level (a.s.l) and depth values.

Figure 13 .
Figure 13.3D diagram of interpolated lignite layer (black color) correlated with lithogical boreholes regarding lithological model with Vertical Exaggeration x4.Vertical axis represent elevation above sea level (a.s.l) values.

Figure 13 .
Figure 13.3D diagram of interpolated lignite layer (black color) correlated with lithogical boreholes regarding lithological model with Vertical Exaggeration x4.Vertical axis represent elevation above sea level (a.s.l) values.

Figure 14 .
Figure 14.3D model of the two hydrolithological units in the targeted area, regarding hydrauli conductivity m/day.Vertical axis represent elevation above sea level (a.s.l) values.

Figure 14 .Figure 15 .
Figure 14.3D model of the two hydrolithological units in the targeted area, regarding hydraulic conductivity m/day.Vertical axis represent elevation above sea level (a.s.l) values.Mining 2022, 2, FOR PEER REVIEW

Figure 15 .
Figure 15.W-E Hydrolithological cross-section on selected boreholes with Vertical Exaggeration x2.Vertical axis represent elevation above sea level (a.s.l) values.

Table 1 .
Lithological units that correspond in the borehole database.

Table 3 .
Hydrolithological units based on the lithological model.