Geological model, deposit block model and reservoir model report

Gautier Laurent, Jean-Jacques Royer, Fabrice Golfier, Lev Filippov, Yulia

Gruzdeva, Maelle Vergnaud, Caroline Izart, Benedicte Lechenard

► To cite this version:

Gautier Laurent, Jean-Jacques Royer, Fabrice Golfier, Lev Filippov, Yulia Gruzdeva, et al. Geological model, deposit block model and reservoir model report. De-170313-0030, Horizon H2020 Projects. 2017, pp.68. insu-05229618

HAL Id: insu-05229618

https://insu.hal.science/insu-05229618v1

Submitted on 29 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
BIO More
An Alternative Mining Concept

Geological model, deposit block model
and reservoir model report

Deliverable no. 2.2 WP2

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>
**Author:** Gautier Laurent

**Organization:** CNRS

**Contributors:**

Jean-Jacques Royer, Fabrice Golfier, Lev Filippov,

Yulia Gruzdeva, Maelle Vergnaud, Caroline Izart, Benedicte Lechenard

**Checked by:**

Name: CNRS

Date: 07/03/2017

Signature:

<figure><img src="image_4.png" /></figure>

**Approved by:**

Name: DMT

Date: 08/03/2017

Signature:

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>
<table><tr><td>Due date of Deliverable</td><td>Project months 15</td></tr><tr><td>Actual Submission Date</td><td>2017-03-13</td></tr><tr><td>Start Date of Project</td><td>2015-05-01</td></tr><tr><td>Duration</td><td>14 months</td></tr><tr><td>Deliverable Lead Contractor</td><td>CNRS</td></tr><tr><td>Revision</td><td>Version 2</td></tr><tr><td>Last Modifications</td><td>2017-03-08</td></tr><tr><td>Nature</td><td>R</td></tr><tr><td>Dissemination level</td><td>PU</td></tr><tr><td>Public Summary enclosed</td><td>No</td></tr><tr><td>Reference / Workpackage</td><td>WP 2</td></tr><tr><td>Digital File Name</td><td>De-170313-0030 - D2.2<br>Geological model, deposit block model and reservoir model report</td></tr><tr><td>Document reference number</td><td>De-170313-0030</td></tr><tr><td>No of pages</td><td>68</td></tr><tr><td>Keywords</td><td>Geological models</td></tr><tr><td>In bibliography, this report should be cited as follows:</td><td>n/a</td></tr></table>

<figure><img src="image_7.png" /></figure>
## Content

<table><tr><td>Content</td><td>5</td></tr><tr><td>List of figures and tables</td><td>7</td></tr><tr><td>Executive summary</td><td>9</td></tr><tr><td>1. Introduction</td><td>11</td></tr><tr><td>1.1. Context of study</td><td>11</td></tr><tr><td>1.2. Selected site and criteria for developing methodology</td><td>12</td></tr><tr><td>1.3. Modelling scales and resolutions</td><td>12</td></tr><tr><td>1.4. Geological context</td><td>13</td></tr><tr><td>1.4.1. Geological Background</td><td>13</td></tr><tr><td>1.4.2. Mineralization</td><td>14</td></tr><tr><td>1.4.3. Tectonic Settings</td><td>15</td></tr><tr><td>2. 3D modelling of geological structures</td><td>17</td></tr><tr><td>2.1. Method for modelling geological structures</td><td>17</td></tr><tr><td>2.2. Available geological datasets</td><td>17</td></tr><tr><td>2.3. Stratigraphic column</td><td>18</td></tr><tr><td>2.4. Regional model: North Sudetic Trough</td><td>19</td></tr><tr><td>2.5. Simplified regional model: North Sudetic Trough</td><td>19</td></tr><tr><td>2.6. Example of district scale modelling</td><td>23</td></tr><tr><td>2.6.1. 3D Modelling</td><td>24</td></tr><tr><td>2.6.2. Fault Network Modelling</td><td>24</td></tr><tr><td>2.6.3. Horizon Modelling</td><td>25</td></tr><tr><td>2.6.4. Honouring well markers and Editing Fault Throw</td><td>26</td></tr><tr><td>2.7. Example of well field scale modelling</td><td>28</td></tr><tr><td>3. Construction of 3D geological grids</td><td>29</td></tr><tr><td>3.1. Geological grid of the simplified regional model</td><td>29</td></tr><tr><td>3.2. Geological grid of the district scale model</td><td>30</td></tr><tr><td>4. Physical property and grade simulation</td><td>32</td></tr><tr><td>4.1. Property modelling strategy</td><td>32</td></tr><tr><td>4.2. Simulation of porosity and permeability</td><td>32</td></tr><tr><td>4.2.1. Estimation of the distribution of porosity property by region</td><td>32</td></tr><tr><td>4.2.2. Estimation of the spatial variability of porosity property</td><td>35</td></tr><tr><td>4.2.3. Estimation of permeability distribution</td><td>41</td></tr><tr><td>4.2.4. Estimation of the spatial variability of permeability</td><td>41</td></tr><tr><td>4.2.5. Sequential Gaussian Simulation of properties</td><td>43</td></tr><tr><td>4.3. Simulation of Copper grade</td><td>48</td></tr><tr><td>5. Hydrogeological model</td><td>51</td></tr></table>

<figure><img src="image_8.png" /></figure>
<table><tr><td>5.1.</td><td>Characterisation and simplification of aquifers</td><td>51</td></tr><tr><td>5.2.</td><td>Boundary conditions</td><td>51</td></tr><tr><td>5.3.</td><td>Preparations for fluid flows modelling</td><td>53</td></tr><tr><td>5.4.</td><td>Export of flow grid on Eclipse for the fluid flows modelling</td><td>54</td></tr><tr><td>6.</td><td>Partner communication and model export</td><td>55</td></tr><tr><td>7.</td><td>Discussion and conclusions</td><td>56</td></tr></table>

<figure><img src="image_9.png" /></figure>
List of figures and tables

<table><tr><td>Figure 1: Main tectonic structures of North-Sudetic area</td><td>12</td></tr><tr><td>Figure 2. Generalized ore zone stratigraphy (Poland)</td><td>14</td></tr><tr><td>Figure 3: Block-diagram of deposition of sediments of Lubin region</td><td>16</td></tr><tr><td>Figure 4: South view of the new available dataset</td><td>18</td></tr><tr><td>Figure 5: Stratigraphic column used in the structural model.</td><td>18</td></tr><tr><td>Figure 6: North Sudetic Trough Model view from South-West</td><td>19</td></tr><tr><td>Figure 7: North Sudetic Trough Model view from top</td><td>19</td></tr><tr><td>Figure 8: Resulting structural model after fault network simplification</td><td>20</td></tr><tr><td>Figure 9: Model with layers from Rotliegendes to Upper Cretaceous</td><td>21</td></tr><tr><td>Figure 10: Model with layers from Rotliegendes to Keuper</td><td>21</td></tr><tr><td>Figure 11: Model with layers from Rotliegendes to Buntsandstein</td><td>21</td></tr><tr><td>Figure 12: Cross-sections through the simplified regional model</td><td>22</td></tr><tr><td>Figure 13: Interpreted fault network</td><td>26</td></tr><tr><td>Figure 14. Final 3D model of the area</td><td>28</td></tr><tr><td>Figure 15: Final stratigraphic grid for the simplified regional model</td><td>30</td></tr><tr><td>Figure 16: Details of the regional stratigraphic grid.</td><td>30</td></tr><tr><td>Figure 17: Resulting geological grid of the study area.</td><td>31</td></tr><tr><td>Figure 18: Well-log porosity study</td><td>33</td></tr><tr><td>Figure 19: Porosity distribution</td><td>34</td></tr><tr><td>Figure 20: District model porosity measurements</td><td>35</td></tr><tr><td>Figure 21: Experimental variogram (region above the known aquifer)</td><td>37</td></tr><tr><td>Figure 22: Experimental variogram (aquifer region)</td><td>38</td></tr><tr><td>Figure 23: Experimental variogram (region under the known aquifer)</td><td>39</td></tr><tr><td>Figure 24: Histograms of porosity for each region of the district model</td><td>40</td></tr><tr><td>Figure 25: Porosity distribution for the regional model</td><td>41</td></tr><tr><td>Figure 26: Log K - porosity relationship</td><td>42</td></tr><tr><td>Figure 27: Estimated permeability (K)</td><td>42</td></tr><tr><td>Figure 28: Proportion of porosity before the simulation (wells data)</td><td>44</td></tr><tr><td>Figure 29: Proportion of porosity after simulation</td><td>44</td></tr><tr><td>Figure 30: Results of Sequential Gaussian Simulation for porosity</td><td>45</td></tr><tr><td>Figure 31: Results of Sequential Gaussian Simulation for permeability</td><td>46</td></tr><tr><td>Figure 32: Simulation of porosity on the simplified regional model</td><td>47</td></tr><tr><td>Figure 33: Copper grade distribution</td><td>48</td></tr><tr><td>Figure 34: Crossplot of Cu versus Ag in the Kupferschiefer deposit</td><td>48</td></tr><tr><td>Figure 35: Variogram for Copper in the Kupferschiefer sandstone</td><td>49</td></tr><tr><td>Figure 36: Simulated Copper grades in the Kupferschiefer sandstone</td><td>50</td></tr></table>

<figure><img src="image_10.png" /></figure>
<table><tr><td>Figure 37: Identification of the various possible aquifers</td><td>52</td></tr><tr><td>Figure 38: Cross section of aquifers and possible water flow</td><td>52</td></tr><tr><td>Figure 39: District model flow grid loaded in Eclipse flowviz</td><td>54</td></tr><tr><td>Figure 40: The final model in .obj format loaded in Meshlab software</td><td>55</td></tr><tr><td>Figure 41: Cross section number 1 without vertical exaggeration</td><td>58</td></tr><tr><td>Figure 42: Radial well field pattern</td><td>60</td></tr><tr><td>Figure 43: Rectangular well field pattern</td><td>60</td></tr><tr><td>Figure 44: Injection pressure propagation</td><td>64</td></tr><tr><td>Figure 45: Visualisation of failure criterion</td><td>65</td></tr><tr><td>Figure 46: Hydraulic fracture propagation</td><td>66</td></tr><tr><td>Figure 47: Hydraulic fracture propagation</td><td>67</td></tr></table>

<table><tr><td>Table 1: Vertical size and resolution of the simplified regional model</td><td>29</td></tr><tr><td>Table 2: Lateral size and resolution for the simplified regional model</td><td>29</td></tr><tr><td>Table 3: Settings used for the grid modelling in vertical direction and volume of units</td><td>31</td></tr><tr><td>Table 4: Definition of the porosity regions for the regional model</td><td>34</td></tr><tr><td>Table 5. Variogram parameters for porosity in the district model</td><td>36</td></tr><tr><td>Table 6. Values of porosity and permeability used for calibration</td><td>42</td></tr><tr><td>Table 7: Data on aquifers of the study picked up in Poland from Tarkowski et al. (2009)</td><td>53</td></tr><tr><td>Table 8: Estimation of NTG for each stratigraphic unit</td><td>53</td></tr><tr><td>Table 9: Summary of characteristics of both radial and vertical well field design</td><td>61</td></tr><tr><td>Table 10: Geomechanical properties of shale used in hydraulic fracturing simulation</td><td>62</td></tr><tr><td>Table 11: Parameters used for hydraulic fracturing simulation</td><td>64</td></tr></table>

<figure><img src="image_11.png" /></figure>
Executive summary

Geological models of subsurface structures and physical properties are essential ingredients for studying physical processes and preparing for in-situ leaching studies. In this report, we present the modelling strategy to build models for representing the geological structures and for simulating fluid flow and petrophysical properties.

Three classes of models have been generated using GOCAD-SKUA® software, with different scale and level of detail:

* A detailed regional scale model, which contains a model of the geological structures (main faults and horizons) in the North Sudetic Trough area. It is based on models produced in the framework of ProMine project.

* A simplified regional scale model, which comprises a model of geological structures, a stratigraphic grid with petrophysical properties and a hydrogeological model for Eclipse fluid flow simulator. It covers the same region as the detailed regional model but the geological structures have been simplified.

* A district scale model, which covers a smaller region delineated within the regional model. It contains a structural model, a stratigraphic grid with grades and fluid flow properties, and a fluid flow simulation model for Eclipse.

* An example of well field scale model containing petrophysical properties for studying fracture propagation processes.

The different models have been exported in different file formats to be easily shared with BIOMOre partners. To this end, file format converters have been developed based on the RINGMesh library. Structural models are presented in the form of sets of geological surfaces represented by triangulated surfaces. Fluid flow models are represented as stratigraphic grids exported in gocad and eclipse formats.

These models are not to be taken for an accurate representation of structural and petrophysical properties. Because of the limited amount of data, especially for controlling petrophysical properties and grades, assumptions have been taken in the modelling processes. The proposed simulation strategy produces properties that globally honour the statistical data but whose value might locally differ from actual values where it is not directly constrained by observations.

This does not lower the interest of such models as the main objective of this study was to define a modelling strategy and produce a set of models for feeding the development of the modelling tool box in WP2. This only means that additional data should be collected for locally constraining petrophysical properties and structures in preparation for pilot plant design studies.

<figure><img src="image_12.png" /></figure>
<figure><img src="image_13.png" /></figure>

**Graphical abstract:** overview of the models and their relationship with other work package tasks.
# 1. Introduction

## 1.1. Context of study

The work package 2 of BIOMOre project aims at developing tools and methods for scaling up laboratory experiments to pilot and industrial scale. This methodology is to be used as a basis for future prefeasibility studies. In addition to the simulation of physical and chemical processes, this task involves characterising structural and petrophysical specificities of the different targeted deposits. Geological structures influence the spatial distribution of rock properties such as fluid flow properties and chemical composition, which will in turn affect the efficiency of the leaching process. The purpose of task 2.2, and its corresponding deliverable D2.2, is to provide a set of numerical models of underground geological settings that can be used for supporting prefeasibility studies.

At the current stage, the project is not yet ready for scaling up to pilot and industrial test. A prototype is currently developed in the galleries of an underground mine in Poland. As the project goes, new criteria are identified, which might invalidate the sites previously selected. Therefore, this study is not to be taken as part of an actual prefeasibility study. The produced models are instead developed to serve as a common basis for the different tasks involved in the development of numerical simulations for studying in-situ leaching. They provide examples of typical geological configurations that can be adapted for the purpose of developing methodology.

As a matter of fact, available data are not sufficient for producing a reliable model over the whole selected area at the required resolution for planning a pilot test. More data would have to be collected to better characterise the local structures and rock properties.

The study presented in this report illustrates four main stages of the modelling process:

1. Building structural models at various scales (Section 2).

2. Populating grids with petrophysical and chemical properties (Section 4).

3. Building hydrogeological models (Section 5).

4. Exporting and sharing the structural models among partners (Section 6).

A short review is first given in Section 1.4 to present the structural evolution of the region as well as the context in which the mineralisation formed.

The various structural models presented in this report have been produced with GOCAD-SKUA<fnref n="1" /> 15 with the help of some research plugins from the Gocad Research Consortium<fnref n="2" />. The coordinate system used in this study is WGS84/UTM33N.

<figure><img src="image_14.png" /></figure>

<footnote n="1">¹ http://www.pdgm.com/products/skua-gocad/</footnote>

<footnote n="2">² http://ring.georessources.univ-lorraine.fr/</footnote>
## 1.2. Selected site and criteria for developing methodology

Deliverable 2.1 has selected the North Sudetic Trough as region of interest for this study (Figure 1). The region of Weisswasser in particular is of interest for building district scale models as additional data have been recently collected (Section 2.1). In this report we have followed the region delineated in D2.1 report. This is mainly because additional data are available at this location, but this is not to be taken as an indication that a pilot plant would or should be implemented there.

## 1.3. Modelling scales and resolutions

Two different scales have been considered in this study (Figure 1):

* Regional scale: covers the whole North Sudetic Trough (approximately 140km by 40 km) for which a detailed model and a simplified model have been produced.

* District scale: scale that would typically be necessary to consider for planning a pilot plant. The chosen district is centred on a recent seismic survey that would help constraining a higher resolution model as compared to the regional one. The district model is about 19km by 12km.

<figure><img src="image_15.png" /><figcaption>Figure 1: Main tectonic structures of North-Sudetic area with designation of copper controlling faults (adapted from Hartsch, 2015).</figcaption></figure>

<figure><img src="image_16.png" /></figure>
## 1.4. Geological context

### 1.4.1. Geological Background

The Kupferschiefer covers an area of approximately 600,000 km² in Central Europe and ranges from England through the Netherlands, Germany to Poland, the Baltic countries and Belarus. It is dated from Upper Permian (Zechstein) (Mistra, 2000). The Kupferschiefer ("copper shale") is a stratigraphic term for a thin layer of finely laminated marine, pyritic, bituminous marl with dispersed sulphidic minerals such as chalcocite, bornite and chalcopyrite (Kutschke et al., 2015).

Figure 2 shows a stratigraphic vertical section of the Lubin region associated with the Kupferschiefer mineralization. The Kupferschiefer (T1) is underlined by the Lower Permian Rotliegendes (red sandstone) that can be characterized by terrestrial silicoclastic sediments consisting of sandstones, silt and mud-stones. The lower part of the Rotliegendes is composed of continental redbeds containing volcanics and in the upper part of sandstones. The uppermost unit of the Rotliegendes sandstones, the Grauliegend or Weissliegend, may reach a thickness of more than 30 m and has been partly eroded during the transgression of the Zechstein Sea. The deposition of the Rotliegendes took place under arid to semi-arid conditions, within the zone of intracontinental rifting at the end of Carboniferous (Vaughan et al. 1989). The block diagram of the sedimentary evolution sequence is given on Figure 3.

In the late Permian, the transgression of the Zechstein Sea from the North and the eustatic rise of the sea level deposited marine shales, dolomites and evaporites while some redbeds were reworked, sometimes bleached due to euxinic conditions of the sea. Therefore, Kupferschiefer is the lowest layer at the base of Zechstein that can be characterized as black shale, deposited randomly on the older rock. The hanging layer, above the Kupferschiefer, is represented by carbonate rocks of the Zechstein Limestone (Ca1) that consists of limestones/dolomites, and corresponds to the Zechstein dolomite and the Zechstein kalk in the German stratigraphic system. The Zechstein Marl (or Mergel in German) (Ca1m) is the clay-rich transitional zone at the basis of limestone. The Zechstein group is divided into several groups showing sequences of evaporation and chemical precipitation. The evaporitic cycles includes Z1-Werra, Z2-Stassfurt, Z3-Leine and Z4-Aller (Zientek et al., 2015).

Lithologic variation within the Kupferschiefer can be explained by paleo-topography of the basin floor: in the deeper parts, it comprises black shale with carbonate-rich layer, while at the edges of basin it is more dolomitic and silty. The thickness of the Kupferschiefer and its marly equivalents change in the north-Sudetic Through from few cm to thickness greater than 2 m (Hartsch 2015). Next period of extension and rifting during the Triassic and Jurassic is probably related with the opening of Tethys Ocean to the south (Vaughan et al., 1989). The deposition of continental marine sediments up to 2,500 m in linear basin took place in the following sequences: Buntsandstein sandstone during the period of the sub-extensional subsidence, then the Muschelkalk and the Keuper marine carbonates and shales within tectonically quieter, thermal subsidence phase (Figure 3).

<figure><img src="image_17.png" /></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 2. Generalized ore zone stratigraphy (Poland) (adapted from Borg et al. 2012).**</figcaption></figure>

At the Jurassic-Cretaceous boundary, the Alpine orogeny that accompanied the closing of the Tethys Ocean caused the inversion and uplift of Lower Silesian basin, creating the Fore-Sudetic block, bounded by NW-SE faults. The uplift is estimated to reach up to 1,000 m in the Lubin region (Poland) that caused the erosion of rocks from Triassic to Cretaceous (Mejia-Herrera 2014). The P-T conditions induce natural hydro-fracturing of the sediments. This naturally increased the rock permeability and activated the circulation of hydrothermal brines. The Permian rocks were tilted 3° to 6° to the northeast, probably during the Alpine orogeny.

Natural hydro-fracturing was an important factor for remobilizing the metals in Kupferschiefer formation, favoured by the replacement of early bacterial pyrite. The deposition of copper in the Weisswasser region is situated in the crossing area of the NNE-SSE and NNE-SSW striking fault systems (Kopp et al., 2012). The tectonic faults of Vascarian and Late Vascarian were reactivated during the Mesozoic Period, giving a possibility to ore bearing brines to flow into the Zechstein sequences. Major source of the metals is considered to be the Lower Permian andesitic and rhyolitic rocks of the North German-Polish Depression. Another source of metal might be earth mantle.

### 1.4.2. Mineralization

The polymetallic Kupferschiefer-ore in the North-Sudetic monocline has a list of various minerals. The major minerals are: chalcocite, bornite, chalcopyrite, covellite, digenite, pyrite, galena, sphalerite, tennantite, native silver, cobaltite, and castaingite. Additionally, minor Cu, Ag, Fe, Ni, Co, Mo, Hg, Bi, Pt and Pd sulphides, arsenides,

<figure><img src="image_19.png" /></figure>

EUROPEAN UNION

This project is funded by the European Union's Horizon 2020 research and innovation programme under grant agreement No 642456.
diarsenides, sulphosalts, thiosulphates, arsenates, noble metals and alloys containing Au, Ag, Hg and Pb are presented in minor quantities (Piestrzyński & Sawlowicz, 1999; Marion et al., 2015).

Piestrzyński and Sawlowicz (1999) determined three stages of mineralization:

1. Syn-sedimentary to early diagenetic bacterial pyrite mineralization phase sustained by the presence of framboidal pyrite;

2. Then replacement phase followed by the formation of the Cu (or Cu-Fe) sulphide related to hydrothermal fluids;

3. The last stage during which Au, other precious metals and significant quantities of Cu have been remobilized.

### 1.4.3. Tectonic Settings

Tectonic conditions are important for understanding the geological evolution of the study area and for building the 3D model. A brilliant summary of the tectonic conditions is provided by Hartsch (2015).

At the North-Sudetic Through area, the faults predominantly strike in NW-SE and NNE-SSW. The shear fault system was formed due to the stress field with E-W maximum extension during eastward drift of European plate to the African plate 300 Ma ago (Borg et al., 2012). The shearing facilitated the extrusion of magma and the extension of volcanic rocks that lead to the deposition of rhyolitic and basaltic to andesitic volcanic rocks as well as some clastic sediments of the Lower Rotliegend.

The Mulkwitz structure is located between Spremberg and Weisswasser, and was developed during the Saxonian movements (Figure 1). It can be described as a horst, developed during the SW-NE compression. The bedding of the mineralized layers at the Zechstein basis are much more tectonically disturbed and more complicated than in other regions of the Central European Copper Belt.

The local uplift of blocks between the Mulkwitz Structure and the Lusatian Block, induced by crustal heating and extension, caused erosion of the pre-Permian rocks. This event reduced the thickness of Upper Rotliegend. The Mesozoic sequences were intensively eroded due to the uplift and Tertiary sediments occur discordantly on the top of the Cretaceous to Triassic formations.

<figure><img src="image_20.png" /></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 3: Block-diagram of deposition of sediments of Lubin region with geological time-scale with interpreted local kinematic regime (adapted from Mejia-Herrera, 2014).**</figcaption></figure>

The variation of the base of the Rotliegend shows no variation in a NW-SE direction in the study area, however, perpendicular to the Main Lusatian Fault, the morphological structure shows lowering of the base of Rotliegend up to 100 m at the SW flank of the Mulkwitz structure and about 150-200 m below NE flank.

The NE-SW fault at the Mulkwitz Structure has been interpreted as downthrown faults by (Hartsch, 2015) with displacement up to 160 m. There are two main regional faults in the NW and SE direction of the structure: the Cretaceous Turkendorf Graben and the Weisswasser fault zone. During the lower Zechstein and from the Cretaceous to the Tertiary, these structures were reactivated.

The crossing of the Weisswasser Fault zone with the SE part of the Mulkwitz Structure has been influenced by movements in different geological periods. As noted by Hartsch (2015), it is possible to see the extreme variation of the thickness of Zechstein or Tertiary sediments over the short distances of 100 m, independently from the regional tends.

<figure><img src="image_22.png" /></figure>
## 2. 3D modelling of geological structures

### 2.1. Method for modelling geological structures

Three-dimensional structural models at various scales have been modelled thanks to SKUA's structural workflow. This tool relies on a powerful implicit modelling approach, which represents faults and stratigraphic units as iso-surfaces of global scalar fields. This method simplifies complex 3D mesh operations required during the structural modelling process and makes it possible to consider several conformable stratigraphic units at once, which improves the quality of data integration. The modelling stages that are specific to each scale and application have been detailed in the corresponding sections.

### 2.2. Available geological datasets

To build the various 3D model of the regional and district area, the following datasets have been used:

* A geological map of the area (in JPEG format) used for identifying new faults that were missing in the previous model.

* Seismic profiles acquired by the Geophysics GGD Company in 2011, comprising four 2D vertical seismic profiles, about 44 km in length and 3 km in depth (in pdf format). The seismic data provides a global architecture of the underground geological structures, more accurate contacts between horizons, and were used for creating lineaments for building the horizons.

* Four drill holes containing geological and geophysical log interpretations and including some geochemical analysis for some of them (in pdf format).

* A 3D model conducted by Mejia-Herrera (2014) as a part of his PhD thesis in the framework of the ProMine project. The model is in coordinate system WGS84/UTM33N.

The 3D model of North-Sudetic Through covers an area of about 700 km². It was created based on 288 wells provided by GEOS and KGHM Kupfer, a geological map for identifying the out-cropping formations, fault traces digitalized from a 2D structural map, interpreted cross-sections, and geological maps showing the contours of the Zechstein that were provided by KGHM Polska Miedz S.A. (Mejia-Herrera, 2014).

The district zone was modelled based on the regional model and upgraded accounting for a higher number of details. This detailed area is located at the intersection of the Mulkwitz Structure and the Weisswasser fault zone. It is oriented in such a way to include the newly acquired available dataset from the seismic survey, the four new drill cores and majority of older drillholes (Figure 4). In the end, the district model incorporates 198 drillholes.

<figure><img src="image_23.png" /></figure>
<figure><img src="image_24.png" /><figcaption>**Figure 4: South view of the new available dataset showing the Volume of Interest of the studied area, the four seismic sections and four wells available from KGHM Kupfer.**</figcaption></figure>

<figure><img src="image_25.png" /></figure>

(Continued on next column)

<figure><img src="image_26.png" /></figure>

(Continued on next column)

<figure><img src="image_27.png" /><figcaption>**Figure 5: Stratigraphic column used in the structural model.**</figcaption></figure>

## 2.3. Stratigraphic column

A complete stratigraphic column (Figure 5) has been defined in GOCAD-SKUA based on the different well markers and horizons available in the datasets. It comprises three major unconformities at the base of the Triassic, upper-cretaceous and Rotliegendes layers. They cause discontinuities in the stratigraphy and are in place associated with the disappearance of eroded horizons. For example, Keuper and Muschelkalk layers are only present in a limited of zones where data support their occurrence. In other places, the interpolation process suggests these layers have been eroded as they are computed to be placed above the corresponding interpolated unconformity.

<figure><img src="image_28.png" /></figure>
## 2.4. Regional model: North Sudetic Trough

<figure><img src="image_29.png" /><figcaption>**Figure 6: North Sudetic Trough Model view from South-West vertical exaggeration x3 (with Gocad-Skua 2009.4)**</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Figure 7: North Sudetic Trough Model view from top (with Gocad-Skua 2009.4)**</figcaption></figure>

A detailed model of the regional North Sudetic Trough has first been produced based on the existing model built in the framework of ProMine (Figure 6 and Figure 7) (Mejia-Herrera, 2014). In practice, little additional work has been needed for this stage. It consisted mainly in checking and editing the structural and stratigraphic information attached to the different parts of the model to ensure its consistency.

The lateral extension of the model is roughly 140km by 40km. In the vertical direction, it extends from +790m to -2200m. It shows a complex framework of faults.

## 2.5. Simplified regional model: North Sudetic Trough

A simplified regional model has been generated for the purpose of running large scale physical simulations for fluid flow or stress state modelling. Simplifications consisted mainly in simplifying the fault network and accordingly updating the stratigraphic horizon interpolation.

In the end, 8 major faults have been kept in the study area. The simplification is done in the “Structure and stratigraphy” SKUA workflow. Fault segments that are close and parallel have been gathered into one large fault. Smaller strike slip faults affecting the geometry of the large ones have been removed and the geometry of the large faults has been made continuous. Shallow faults or faults close to the boundary of the model have been removed as they would have little impact on simulations.

<figure><img src="image_31.png" /></figure>
<figure><img src="image_32.png" /><figcaption>**Figure 8: Resulting structural model after fault network simplification.**</figcaption></figure>

The well data have been inspected for quality control. We checked that the resulting fault network was consistent with the location of horizons. Some faults have been extended in the step "Fault Contact Modelling", Sued3 and FailleN2, so they touch the Hoyerswerda fault and respectively the Weisswasser fault. Appropriate structural features have been assigned to each segment of fault. In practice, each segment, even belonging to the same regional structure, needs to be described by its own geological feature. This is one of the limitations of the software.

We wanted to be more specific about the extension of faults at depth. Some faults appear on the geological map so we suppose that they cross all layers of the model. For faults which appear partially on the geological map, we supposed they cross all layers to the limit Cretaceous/Tertiary. Faults which don't appear on the geological map were supposed to cross all layers to the limit Keuper/Cretaceous.

In the end, the model consists of 8 faults. They account for the principal directions of the main faults of the initial model and define similar regional blocks. Geological horizons were updated to take account of the simplified fault model. Any gap or throw in the original horizon were replaced by a continuous piece of surface by interpolation. The resulting model is water-tight by construction. The resulting model is illustrated in Figure 8 to Figure 12 and appendix 1.

<figure><img src="image_33.png" /></figure>
<figure><img src="image_34.png" /><figcaption>Figure 9: Model with layers from Rotliegendes to Upper Cretaceous</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Figure 10: Model with layers from Rotliegendes to Keuper</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Figure 11: Model with layers from Rotliegendes to Buntsandstein</figcaption></figure>

<figure><img src="image_37.png" /></figure>
<figure><img src="image_38.png" /><figcaption>a. Location of cross sections</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>b. Cross section AA' oriented WNW-ESE</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>c. Cross section BB' oriented SW-NE</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>d. Cross section CC' oriented SW-NE</figcaption></figure>

**Figure 12: Cross-sections through the simplified regional model.**

<figure><img src="image_42.png" /></figure>
## 2.6. Example of district scale modelling

The geological map and the seismic profiles were imported as 2D Voxets. The seismic traces sections were created in Gocad/SKUA from the coordinates of the measuring stations recorded in column based file. Subsequently, the Voxets were fit and aligned along the seismic traces using the Gocad GeoToolbox Plugin (Ring.georessources 2016). The geological map was georeferenced using three points based on the well locations.

The four new wells referred as WW01, WW02, WW03, WW04 were created from the log data with markers attached to the bottom of the corresponding stratigraphic units. They assumed to be drilled from the surface; therefore, the elevation coordinates were derived from the Digital Elevation Model DEM of this area. The wells are assumed to be vertical due to the fact that available data showed a very limited deviation from the vertical direction.

The chemical assays, comprising Cu, Pb, Zn and Ag and some other elements, including REE trace and major elements Au, Pt, Pd, Mo, Ni, Co, Mn, As, Fe, U, Th, Sr, Cd, Sb, Bi, V, Ca, P, Cr, Mg, Ba, Ti, Al, K, W, Zr, Ce, Sn, Y, Nb, Ta, Be, Sc, Li, S, Rb, Cd, Re, Mn, Ga, Cs, Ge (in ppm or %), were collected in an Excel™ sheet, and imported as a column-based file for Cu, Pb, Zn and Ag. A surface copper density (in kgCu.km⁻²) was then calculated for each drill core assuming two cut-offs at 0.1 and 0.3 % in Cu, respectively, for delineating the mineralized zone. Based on these four drill holes, the samples of the Weisswasser area exhibit in general a lower concentration in metal content compared to the samples from the Polish mines for which the average values being in the range of 1.4% for Cu (Zientek et al., 2015). High content of calcium observed in the mineralized zone which may deteriorate the efficiency of the leaching process.

In addition, the stratigraphic column used for the modelling in the previous study by Mejia-Herrera (2014) was modified to account for details observed on the new four wells of the Weisswasser area. The geological cross-sections of the Spremberg area come from a nearby study conducted by Kopp et al. (2006) on the NE of Weisswasser showing evidences of unconformable contacts, in particular due to erosion and baselap, for the Strassfurt and the Upper Cretaceous units, respectively. Moreover, drill holes performed in the volume of interest (VOI) of the study area do not show evidences of the presence of the Keuper formation; therefore, this horizon which is present in the stratigraphic column was omitted in the model (Figure 5).

The horizon traces were drawn from the vertical seismic profiles. They include the base of the Kreide, The top Zechstein and the Base of the Zechstein horizons as drawn by Geophysic GGD. They correspond to the Low_Creataceous, Z4_Allers, Z1_Zechsteinkonglomerat, respectively. These lines were assigned to the horizons to be used in the workflow. It should be mentioned that the seismic data is less precise than the well data. In case of inconsistency on the depths with the well markers values, seismic data was ignored. Generally, the depths picked on the seismic vertical sections showed good fit with the well markers. Only, one drill hole (WW1102) located in the vicinity of the seismic profile WW1103_seismic shows some inconsistency.

<figure><img src="image_43.png" /></figure>
A quick-look at the newly acquired dataset has revealed some inconsistencies in the previous 3D model that justify a new reinterpretation. For example, strong variations in the markers elevation of the same unit suggest the existence of faults that were not accounted for in the previous model. This is why a detailed analysis of the fault system was carried out. This analysis suggested the presence of new faults as shown on the geological map. These new faults were introduced in the new 3D model. Since the dipping angles of the faulting are unknown, they were assumed to be vertical and of the same vertical and lateral extend as the-neighboured faults.

2.6.1. 3D Modelling

The 3D model is built using the *Structure and Stratigraphy Workflow* available in Gocad/SKUA. The workflow provides an implicit modelling approach with a geological grid as an output. The workflow is organized following the subsequent steps: (i) data preparation and addition, (ii) fault network building, (iii) horizon building, (iv) honouring well data, and (v) building the final geological grid. The advantages of this modelling approach are that the model can be dynamically upgraded when the newly acquired data are available or for accounting for new structural interpretations.

A data import is the preliminary step of the workflow; it includes a conceptual data interpretation, such as a choice for the fault type (normal, reverse or strike-slip), the presence or not of unconformities, of eroded faults, etc. Therefore, the geo-chronological workflow implies that the user have to investigate in advance all the information necessary to achieve these steps. However, the lack of information on complex data such as faults is often a limiting factor for the interpretation. Working in 3D space without geological cross-sections, and/or 3D seismic cubes makes difficult to understand the spatial trends of the subsurface structures. In this particular case, considering the faults network, no information was available on the fault throw types. Thus, it was decided not to specify the fault throw, except in the obvious case of listric faults or when it was possible to identify reverse faults.

The Volume of Interest (VOI) was chosen in order to include the newly acquired seismic cross-sections and the well data and to fit and aligned the NE-SW direction of the main fault network. The VOI was built using a rectangular frame with the top surface set at the elevation of 150 m and the bottom surface set at -2,200 m, with extends of 18,965 m in length and 11,845 m in width. Existing faults and horizons outside of VOI from the previous model were eliminated (Figure 4).

2.6.2. Fault Network Modelling

The fault network was upgraded as newly acquired data showed that the previous structural interpretation was not fully correct. The goal of this step is to generate the fault blocks on which the tetrahedral mesh is propagated, respecting the contact relationships. The workflow automatically checks the validity of contacts and points out mistakes. Therefore, the building of the faults network is time consuming as it includes manual editing of the contacts between faults, faults shapes, etc. The stratigraphic unit are assumed to be continuous within the fault blocks that is used for horizon building on the next step.

To do this efficiently, the geological context of the studied area have to be fully understood. By looking on the geological map only it is challenging to determine the

<figure><img src="image_44.png" /></figure>
relationships between faults (fault order). Generally, the area is characterized by NE-SW and SSE-NNW trending faults. According to the tectonical evolution discussed in the literature survey, the oldest identifiable structures are NE-SW trending faults with steep dipping, caused by a SE-NW compression during the Saxonian orogeny. Therefore, the NE-SW structures were formed as first order faults and have been cut by another set of faults in the SSE-NNW direction. The previous model is in agreement with this interpretation.

New faults have been included in the model. The motivation for adding new faults was the mismatch between the geological markers of the same kind observed at the new well location; the location of new faults was derived from the geological map. Five new surfaces were created with prefix ‘d’ in their name since they were shown as dashed faults on the geological map.

In the central part of the study area, there is the Weisswasser fault zone bounded on the East side by the fault 221 and from the West side by the Weisswasser faults, which are interpreted as normal faults with a SE-NW orientation and vertical dip. Despite the fact that the area is interpreted as a through zone, the fault down throws are not significant. The area is strongly faulted in the lateral direction. The maximum vertical displacement about 200 to 300 m was observed on the 221 fault (Hartsch, 2015) indicating a down-throw at the eastern border of the study area. The initial and current fault networks are shown in Figure 13

2.6.3. Horizon Modelling

The horizons building algorithm follows the sequence of the stratigraphic column from top to bottom, using assign data to the horizon as control points, and simultaneously obeying the faults contact relation for all horizons. The workflow allows the selection of different horizons in previewed mode that helps to indentify which parameters can cause inconsistencies in the model. Gradually adding the underneath horizons and adding/ignoring initial data (for example, parts of the seismic lines or additional lines) helps in finding solutions to inconsistencies. When all the horizons are geologically consistent, all horizons are built in one step.

Horizon modelling of the study area was extremely tedious. This was due to: (i) the
dense faulting, (ii) to significant difference in markers elevation between different parts
of the model, and (iii) to relatively small thickness of some units compared with lateral
extend horizon thickness, in particular, for the oldest Zechstein formation. Some
problems were related to the absence of data in the some of the blocks. To proceed
with the workflow and to fix these flaws, some additional lines for better defining the
horizons were added manually in some locations. They were used further as
constrains to facilitate the horizon modelling. For example, in the SE of the Mulkwitz
structure, there is a significant down through observed at the drill hole CuSp_28_60
that was considered at the beginning by the software as an outlier. The seismic data
is located in the opposite part of VOI, and could not help the interpretation in this case.
Therefore, this block was modelled including manually control lines consistent with the
observed down through.

<figure><img src="image_45.png" /></figure>
<figure><img src="image_46.png" /><figcaption>**Figure 13: Interpreted fault network initial (top) (Mejia-Herrera, 2014) and current model (bottom).**</figcaption></figure>

Another recurrent problem is the high variability observed on the thickness of some irregular horizons (called uneven horizons in the Gocad/SKUA Workflow). This implies that layers may have different sediment rates resulting in different thickness. Since the whole horizons pile presents slightly dipping set of sediments on the flank of Mulkwitz structure, it is acceptable to have slight variation in thicknesses. However, when setting this constrain, the interpolation algorithm does not precisely fit the hard control data (markers); consequently, the next step of the workflow "honouring markers" is not performed properly. For the current model, the uneven thickness constraint was set on the Q (Quaternary), Tr (Triassic), and Strassfurt units, not specifying the direction of the grid (azimuth) in the horizons.

During the modelling step, numerous unconformities have been identified: Low and Upper Cretaceous and Muschelkalk are partially eroded; Buntsandstein, Aller and Leine are deposited in baselap mode (overlying the Strassfurt horizons). Concerning the Zechstein series, their lineaments are quite similar all over the study area. Differences in thicknesses exist in the NE and SW parts of the study area: the NE part contains approximately 380 m thick of sediment, and the SW part 270 m.

#### 2.6.4. Honouring well markers and Editing Fault Throw

<figure><img src="image_47.png" /></figure>
As control lines picked on the vertical seismic sections (or additional control lines introduced manually in some locations) are considered as control lines, some inconsistencies may occurred with the well markers which may not match perfectly the horizons. This model will be used for quantifying the resources in place and to drive fracturing simulation in some specific mineralized layers. Therefore, it is important to respect as much as possible the position of horizons observed on the exploration wells in order to start with an *as-much-as possible* "true" 3D model. So an additional manual step was performed in order to force the horizon surfaces to pass accurately though the well markers.

However, this manual fitting step results in creating some geologically unrealistic rough disturbances in the horizon surfaces, especially nearby wells and seismic cross-sections that have been propagated to underneath horizons. The improvements of the horizons surfaces require a lot of manual editing in the *editing fault throw* step of the Workflow. It allows modifying the surfaces by redrawing manually the horizon displacement lines on each side of the faults. The main goal of this step is to avoid artificial occurring of fault throw, i.e. the fault throw should be constant, and have a constant displacement magnitude within a given layer.

Redrawing a fault contact lines can lead to marker mismatch, especially for boreholes located close to the editing fault. To overcome this problem, a fitted through markers surfaces was used to generate additional control lines. Then, the *horizon building* step was reiterated turning on active the *maximum data fit* constrain, and turning off the *marker fit* constrain. The mismatch is still present, but significantly reduced. Accounting for the fact that some geological units are very thin, reaching values as low as 20 cm for the Kupferschiefer shales, these mismatches can be an acceptable approximation. The final 3D model of the underground structures of the Weisswasser area is shown in Figure 14.

Another way to build a reliable model for fracture simulation with little manual editing work is to reduce the size of the model. For example, only relevant horizons could be included in the VOI to create a final geologic grid.

<figure><img src="image_48.png" /></figure>
<figure><img src="image_49.png" /><figcaption>**Figure 14.** Final 3D model of the area built using the Structure and Stratigraphy workflow, representing the faults network with all horizons. Vertical exaggeration × 2.</figcaption></figure>

## 2.7. Example of well field scale modelling

After a model has been generated for the scale of a particular targeted district (Section 2.6), more accurate models would have to be constructed for better constraining the implantation of a possible pilot plant. This requires models at a finer scale and higher resolution. Appendix 2 provides a preliminary study for such a well field scale model.

<figure><img src="image_50.png" /></figure>
## 3. Construction of 3D geological grids

Three-dimensional geological grids have been created for the simplified regional model (Section 3.1) and the district model (Section 3.2). The geological grids are built fitting the fault network structure and the horizon architecture. The grid provides the geostatistical framework for further subsequent petrophysical modelling. They were generated with the *structural grid building* workflow of Gocad/SKUA. Appropriate vertical resolution has to be chosen for each stratigraphic unit.

### 3.1. Geological grid of the simplified regional model

This grid has been built considering the whole stratigraphic pile from the quaternary sediment cover to the Rotliegend formation. The dimensions and resolution of the grid are summarized in the Table 1 and Table 2. The total number of cells is 208 times 79 times 239, hence 3'927'248 cells. The resulting model is illustrated in Figure 15 and Figure 16.

<table><caption>Table 1: Vertical size and resolution of the simplified regional model.</caption><thead><tr><td>Stratigraphic Unit</td><td>Volume (m³)(.10⁹)</td><td>Number of cells</td><td>Cell Thickness</td></tr></thead><tbody><tr><td>Quaternary</td><td>357.959</td><td>12</td><td>19.94</td></tr><tr><td>Tertiary</td><td>212.122</td><td>10</td><td>20.975</td></tr><tr><td>Upper Cretaceous</td><td>689.889</td><td>58</td><td>20.0.35</td></tr><tr><td>Lower Cretaceous</td><td>206.228</td><td>37</td><td>20.148</td></tr><tr><td>Keuper</td><td>108.511</td><td>54</td><td>19.844</td></tr><tr><td>Muschelkalk</td><td>249.539</td><td>8</td><td>19.884</td></tr><tr><td>Buntsandstein</td><td>1,701.505</td><td>24</td><td>19.884</td></tr><tr><td>Aller</td><td>646.1</td><td>8</td><td>20.058</td></tr><tr><td>Leine</td><td>645.974</td><td>8</td><td>19.659</td></tr><tr><td>Strassfurt</td><td>86.782</td><td>1</td><td>20.5</td></tr><tr><td>Oberer Werra Anhydrit</td><td>82.177</td><td>1</td><td>19.105</td></tr><tr><td>Aeltestes Steinsalz</td><td>205.107</td><td>2</td><td>23.371</td></tr><tr><td>Unterer Werra Anhydrit</td><td>495.423</td><td>5</td><td>21.139</td></tr><tr><td>Zechsteindolomit</td><td>273.38</td><td>3</td><td>18.702</td></tr><tr><td>Zechsteinkalk</td><td>17.003</td><td>1</td><td>3.438</td></tr><tr><td>Zechsteinmergel</td><td>8.336</td><td>1</td><td>1.691</td></tr><tr><td>Kupferschiefer</td><td>3.136</td><td>1</td><td>0.635</td></tr><tr><td>Zechsteinconglomerat</td><td>271.523</td><td>3</td><td>17.803</td></tr><tr><td>Grauliegendes</td><td>66.564</td><td>1</td><td>13.061</td></tr><tr><td>Rotliegendes</td><td>59.07</td><td>1</td><td>11.548</td></tr></tbody></table>

<table><caption>Table 2: Lateral size and resolution for the simplified regional model.</caption><thead><tr><td></td><td>Number of cells</td><td>Cell size (m)</td></tr></thead><tbody><tr><td>I Axis</td><td>208</td><td>700</td></tr><tr><td>J Axis</td><td>79</td><td>600</td></tr></tbody></table>

<figure><img src="image_51.png" /></figure>
<figure><img src="image_52.png" /></figure>

<figure><img src="image_53.png" /><figcaption>Figure 15: Final stratigraphic grid for the simplified regional model.</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>Figure 16: Details of the regional stratigraphic grid.</figcaption></figure>

## 3.2. Geological grid of the district scale model

<figure><img src="image_55.png" /></figure>
The grid was built only on the five following facies including the dolomites, carbonates (kalk), marl (mergel), shales (kupferschiefer) and sandstones (conglomerat) (Table 3) because it is designed for resources estimation purposes and for implementing further bioleaching and fracturing simulation which will be focused only on the mineralized zone. Therefore, five regions were created as stratigraphic units in the geological grid.

For the modelling with sufficient resolution, the vertical number of cells in the mineralized zone was increased, especially for the Kupferschiefer horizon. The horizontal resolution is 110 × 100 cells that correspond to about 200 m. The vertical layering of the grid was set to 36 micro layers (vertical cells) unequally spread between the geological units (Table 3). The gross volume of the geological grid (110×100×36) is 25.5 km³ with an average thickness of 100 m. The resulting geological grid is shown in Figure 17.

<table><caption>Table 3: Settings used for the grid modelling in vertical direction and volume of units.</caption><thead><tr><td>Stratigraphic Unit</td><td>Number of vertical cells</td><td>Cell thickness (in m)</td><td>Region Volume (in km³)</td></tr></thead><tbody><tr><td>Zechsteindolomit</td><td>14</td><td>5</td><td>16.3</td></tr><tr><td>Zechsteinkalk</td><td>12</td><td>2</td><td>5.72</td></tr><tr><td>Zechsteinmergel</td><td>6</td><td>2</td><td>2.97</td></tr><tr><td>Kupferschiefer</td><td>3</td><td>0.2</td><td>0.156</td></tr><tr><td>Zechsteinkonglomerat</td><td>1</td><td>1.3</td><td>2.98</td></tr></tbody></table>

<figure><img src="image_56.png" /><figcaption>Figure 17: Resulting geological grid of the study area. Vertical magnification × 2.</figcaption></figure>

<figure><img src="image_57.png" /></figure>
# 4. Physical property and grade simulation

## 4.1. Property modelling strategy

A geostatistical approach has been used to model the physical properties of the different geological layers in the simplified regional model and the district model. For each model, fluid flow properties (porosity and permeability) have been computed. This part describes: (i) the interpretation one can derived from the geophysical data gathered on the boreholes; (ii) the methodology used to populate the geological grid with petrophysical properties (porosity and permeability) using the Reservoir Properties workflow. The main steps can be summarized as follows:

1. Estimation of the porosity ($\phi$) and permeability ($\log K$ and $K$) from petrophysical logs and importing them in the project;

2. Understanding the spatial variability of the petrophysical properties within each geological unit, estimation of the histogram and variogram;

3. Using Sequential Gaussian Simulation (SGS) with ordinary kriging to estimate the porosity distribution in the grid for each geological unit;

4. Using SGS with co-kriging to simulate the permeability distribution conditioned by the porosity, and assuming a log-linear relationship between them.

It should be noted, that while porosity can be directly estimated from the petrophysical logs, the permeability is not so easy to be estimated as there is no direct technique to estimate it from log data. Hence, the permeability is estimated from the log linear porosity - permeability relationship and is calibrated on the experimental data.

Permeability was measured by the UMR7359 GeoRessources laboratory on samples of each main rock types (sandstone, black shale, and dolomite) of the Rudna mines. Experimental results show quite a high dispersion in values for shales. Thus, the permeability values were also fixed from the common range values suggested in the literature. Anisotropy of permeability depends on many factors and is hard to predict, therefore, permeability was assumed to be isotropic.

## 4.2. Simulation of porosity and permeability

### 4.2.1. Estimation of the distribution of porosity property by region

The porosity was estimated from four petrophysical logs (WW1101, WW1102, WW1103, WW1104) by an average point approach. Porosity record is available in all stratigraphic units, but the logging had stopped downward the Zechsteinonglomerat unit. The porosity measurements were estimated from the available diagraphies. For the district model, the values have been directly extracted at a given resolution from the porosities provided with the data. For the regional model, a new estimation has been proposed (Figure 18) and sampled at the scale of the stratigraphic grid.

<figure><img src="image_58.png" /></figure>
<figure><img src="image_59.png" /><figcaption>**Figure 18: Well-log porosity study**
(a) WW01, (b) WW02, (c) WW03, (d) WW04</figcaption></figure>
To extract porosity values in the district model, a sampling step of 10m was used for all units, except for shales and sandstones, where it was reduced to 1m due to the relatively small thickness of those seams. Next, the data was imported as column-based file into the project.

The histogram of the extracted and computed porosity from the four wells is shown in Figure 25 for the simplified regional model, and Figure 20 for the smaller region covered by the district model.

<figure><img src="image_60.png" /><figcaption>Figure 19: Porosity distribution in the different regions of the simplified regional model.</figcaption></figure>

<table><caption>Table 4: Definition of the porosity regions for the regional model.</caption><tr><td>Region 1</td><td>Quaternary</td></tr><tr><td>Region 2</td><td>Tertiary + Upper Cretaceous</td></tr><tr><td>Region 3</td><td>Lower Cretaceous</td></tr><tr><td>Region 4</td><td>Keuper</td></tr><tr><td>Region 5</td><td>Muschelkalk</td></tr><tr><td>Region 6</td><td>Buntsandstein</td></tr><tr><td>Region 7</td><td>Aller</td></tr><tr><td>Region 8</td><td>Leine</td></tr><tr><td>Region 9</td><td>Strassfurt + Oberer Werra Anhydrit</td></tr><tr><td>Region 10</td><td>Aeltestes Steinsalz</td></tr><tr><td>Region 11</td><td>Unterer Werra Anhydrit</td></tr><tr><td>Region 12</td><td>Zechstein + Kupferschiefer + Permian</td></tr></table>

<figure><img src="image_61.png" /></figure>
<figure><img src="image_62.png" /><figcaption>**Figure 20: District model porosity measurements from four wells within the whole geological grid.**</figcaption></figure>

In the district scale model, the histogram of porosity shows distinctly several populations corresponding to rock types. Therefore, the porosity simulation was performed separately within the five units. The porosity values for the different units are summarized in Figure 24 and Figure 25. Low-permeable shales having porosity reaching zero, are overlaid by dolomites, those porosities vary between 0 and 15%. Porosity measurements on underlying sandstones are typically less than 4%, with a peak at 9% observed on the WW04 well. These values are unusually low for sandstone, indicating a possible carbonate digenesis phase that would have plugged the initial porosity confirmed by the petrological study (Marion et al., 2015).

### 4.2.2. Estimation of the spatial variability of porosity property

Geostatistical simulation processes are based on the characterisation of the spatial variability of the quantity being simulated. This is estimated by computing experimental variogram that measure the amount of variability between any two data point with respect to the distance between them. Variograms have been computed in the vertical and horizontal directions for different regions (Figure 19 and Table 4). In the regional model, three regions have been separated: the region above the known aquifer (region 1 and 2, Figure 21), the region of known aquifer (region 3 to 6, Figure 22), and the region under the aquifer (region 7 to 12, Figure 23).

For the district model, the vertical variogram of the porosity was built for all facies. However, all units except dolomites have very few data to build a reliable variogram because of the small number of samples. The parameters for variogram were derived based on the following assumptions:

* The sill of the variogram is assumed to be equal to the variance ($\sigma^2$) taken from histograms;

* To estimate the range, cell size of the grid in the vertical direction was taken into account. Since the UVW transformation was used, the *range is represented*

<figure><img src="image_63.png" /></figure>
as a fraction of the volumetric space. The average thickness of the grid is 100 m;

* A Gaussian type of variogram was chosen, since it provides a 'smoother' distribution due to the up-scaling of the porosity measured each cm into cells of 1 to 5m in the vertical direction. In fact, the Gaussian model has a parabolic rather than a linear behaviour at the origin, which implies slower change in the variability (smaller spatial variability than the spherical model);

The fitted parameters of the variograms are summarized in Table 5.

<table><caption>Table 5. Variogram parameters for porosity in the district model.</caption><thead><tr><th></th><th>Zechstein dolomite</th><th>Zechstein kalk</th><th>Zechstein mergel</th><th>Kupfer schiefer</th><th>Zechstein konglomerat</th></tr></thead><tbody><tr><th>Min value (in %)</th><td>0.04</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td></tr><tr><th>Max value (in %)</th><td>0.15</td><td>0.04</td><td>0.04</td><td>0.03</td><td>0.04</td></tr><tr><th>Sill (in %)<sup>2</sup></th><td>0.0015</td><td>0.002</td><td>0.003</td><td>0.00014</td><td>0.00005</td></tr><tr><th>Range (in %)</th><td>0.150</td><td>0.025</td><td>0.025</td><td>0.005</td><td>0.008</td></tr><tr><th>Type</th><td>Gaussian</td><td>Gaussian</td><td>Gaussian</td><td>Gaussian</td><td>Gaussian</td></tr></tbody></table>

Accurate modelling of the permeability in the reservoir is difficult due to the lack of data. Permeability is one of the parameters that could not be directly estimated from core data. Generally, the permeability derived through a porosity-permeability relationship by building the regression model of the log of permeability (log K) against porosity (φ) for each stratigraphic unit. Afterwards, the permeability is back-calculated from derived log K. The disadvantages of this approach are the variability and the uncertainty in the porosity-permeability scatter plot as a result of sparse sampling leading to a lost in some spatial variability (Boisvert et al., 2012). The resulting data will be used to populate the grid using a Sequential Gaussian Simulation (SGS).

The permeability can be inferred with the following strategy:

1. Estimate Log K from the Log K-φ relationship;

2. Populate Log K in the grid using co-kriging with porosity

3. Compute K from log K applying a property script.

The correlation coefficient between the logarithm of the permeability and the porosity is taken from a bivariate analysis between these properties for each unit.

<figure><img src="image_64.png" /></figure>
<figure><img src="image_65.png" /><figcaption>**Figure 21: Experimental variogram (region above the known aquifer)**</figcaption></figure>

<figure><img src="image_66.png" /></figure>
<figure><img src="image_67.png" /></figure>

<figure><img src="image_68.png" /><figcaption>Figure 22: Experimental variogram (aquifer region)</figcaption></figure>

<figure><img src="image_69.png" /></figure>
<figure><img src="image_70.png" /><figcaption>**Figure 23: Experimental variogram (region under the known aquifer)**</figcaption></figure>

<figure><img src="image_71.png" /></figure>
<figure><img src="image_72.png" /></figure>

**Zechsteindolomit**

Min : 0 %
Max : 15 %
Mean : 3.4 %
St. Dev.: 0.11 %

<figure><img src="image_73.png" /></figure>

**Zechsteinkalk**

Min : 0 %
Max : 12 %
Mean : 5.3 %
St. Dev.: 4.2 %

<figure><img src="image_74.png" /></figure>

**Zechsteinmergel**

Min : 0 %
Max : 11 %
Mean : 4.3 %
St. Dev.: 4.1 %

<figure><img src="image_75.png" /></figure>

**Kupferschiefer**

Min : 0 %
Max : 3 %
Mean : 1.3 %
St. Dev.: 0.57 %

<figure><img src="image_76.png" /></figure>

**Zechsteinkonglomerat**

Min : 1 %
Max : 2 %
Mean : 1.3 %
St. Dev.: 0.57 %

**Figure 24: Histograms of porosity for each region of the district model.**

<figure><img src="image_77.png" /></figure>
<figure><img src="image_78.png" /><figcaption>**Figure 25:** Porosity distribution for the regional model from the interpretation of the four available wells.</figcaption></figure>

#### 4.2.3. Estimation of permeability distribution

The permeability distributions were calculated for the three rock types including dolomite (including Zechsteindolomite, Zechsteinkalk and Zechsteinmergel), shale (Kupferschiefer) and sandstone (Zechsteinkonglomerat). This simplification is done since the measurements for the calibration is provided only for these three rocks. The regression curves for each geological unit were built using typical ranges recommended in the literature for previous similar deposits (Mavko, 2010).

The previous work of the literature suggests a decrease in permeability between geologic facies according to the following trend: sandstone → dolomite → shales. The curves were roughly calibrated with measurements obtained from laboratory experiments. It should be noted that experimental values for shales porosity vary significantly from the observed one (3% against 20%). Therefore, the permeability of shale was reduced to be more consistent with field data. The final values used for the calibration for the porosity and the permeability are summarized in Table 6. The fitted regression curves for the permeability are shown in Figure 26. The estimated Log K and K parameters were imported on new properties on wells in the project from a column-based file. The estimated permeability distribution (in mD) is shown on Figure 27.

In this study we have considered a zonal anisotropy of permeability with a ratio of 1 to 10, i.e. the simulated value of permeability is used directly for PermX and PermY, but is divided by ten for PermZ.

#### 4.2.4. Estimation of the spatial variability of permeability

The variogram is calculated on the porosity and assumed to be identical for the logarithm of the permeability. There is too few data to correctly calibrate this information and the estimation made with the porosity will already give the general trend for the spatial variations.

<figure><img src="image_79.png" /></figure>
<figure><img src="image_80.png" /><figcaption>**Figure 26: Log K - porosity relationship built to estimate the permeability from the porosity for sandstone, dolomite and shales.**</figcaption></figure>

<figure><img src="image_81.png" /><figcaption>**Figure 27: Estimated permeability (K) within the district geological grid (samples number n = 68).**</figcaption></figure>

<table><caption>Table 6. Values of porosity and permeability used for calibration</caption><thead><tr><th></th><th>Sandstone<sup>1</sup></th><th>Dolomite<sup>1</sup></th><th>Shale<sup>2</sup></th></tr></thead><tbody><tr><td>Porosity (in %)</td><td>12.5</td><td>20.0</td><td>1.15</td></tr><tr><td>Permeability (in mD)</td><td>3.6</td><td>1.1</td><td>0.2</td></tr></tbody></table>

<figure><img src="image_82.png" /></figure>

<footnote n="1">¹measured by GeoRessources (Nussipkozhayev, 2015; Marion et al., 2015)</footnote>

<footnote n="2">²taken from literature (Mavko 2010)</footnote>
4.2.5. Sequential Gaussian Simulation of properties

The petrophysical properties were interpolated in the geological grid using the
Reservoir Properties workflow available in Gocad/SKUA. Three continuous properties
were created in the workflow: Porosity, Log K and K. The workflow provides various
geostatistical tools to interpolate the properties distribution. Standard geostatistical
methods for populating a grid with properties include the kriging, the co-kriging and the
Sequential Gaussian Simulation (SGS). The SGS method simulates some realizations
conditioned by well data, and reproducing the histograms and the variograms. Porosity
was simulated by SGS using ordinary kriging, while permeability was simulation by co-
kriging SGS with porosity as secondary variables. SGS allows reducing the smoothing
effect of the kriging and better represents the natural variability of the property (Rossi
& Deutsch, 2013).

The drawback of SGS is that the resulting porosity or permeability field will show less connectivity than the input data. The SGS algorithm steps are:

* To generate a random path within the grid;

* In each cell of the grid, use kriging to estimate a mean and standard deviation based on the surrounding data values;

* Pick at random a value from the corresponding normal distribution and assign the value to that cell;

* Visit each successive node following the random path and repeat the process, taking into account previously simulated values as conditioning data in the kriging process.

Therefore, the SGS requires the knowledge of both the variogram and the histogram
in the stationary domain of investigation. Variogram provides information on the spatial
variability of the studied property against distance between data points, and have to
be modelled for each geologic unit. Modelling variogram requests an estimation of the
ranges and the sill, which turn to be critical input parameters for SGS simulation. The
domain of variation (min and max values) is taken from histograms calculated on the
samples of the corresponding units.

The maximum values used for simulating the porosity was reduced for the
Zechsteinkalk, and Zechsteinmergel due to the fact that the porosity was vertically
averaged on block of 10 m, and the porosity was overestimated. Unknown values for
the ranges were estimated assuming a unit cell. It should be noted that the estimation
from the simulation might not fully reflect the heterogeneity of the area, since very
limited data are available.

The results of the SGS for porosity and permeability are shown on Figure 30, Figure 31 and Figure 32. As expected, each realization shows distinct changes in porosity between layers due to the change in rock types, that may be used for visualization purposes. There is still large degree of uncertainty as a consequence of limited data. As a first step, this is a first approximation of the reality which can be refined as soon as new measurements (wells) will be available.

<figure><img src="image_83.png" /></figure>
After the simulation, the simulated distribution for any region (Figure 28) should correspond to the distribution in the dataset (Figure 29). This is illustrated on the next figures for region 6 of the regional model.

<figure><img src="image_84.png" /><figcaption>Figure 28: Proportion of porosity before the simulation (wells data)</figcaption></figure>

<figure><img src="image_85.png" /><figcaption>Figure 29: Proportion of porosity after simulation on the region 6 of the grid</figcaption></figure>

<figure><img src="image_86.png" /></figure>
<figure><img src="image_87.png" /><figcaption>Figure 30: Results of Sequential Gaussian Simulation for porosity for Zechsteindolomit (A), Zechsteinkalk (B), Zechsteinmergel (C), Kupferschiefer (D), and Zechsteinkonglomerat (E).</figcaption></figure>
<figure><img src="image_88.png" /><figcaption>Figure 31: Results of Sequential Gaussian Simulation for permeability estimated from the Log K for dolomite (A), shale (B) and sandstone (C).</figcaption></figure>
<figure><img src="image_89.png" /><figcaption>**Figure 32: Simulation of porosity on the simplified regional model**</figcaption></figure>
4.3. Simulation of Copper grade

The same simulation strategy has been applied for Copper grades. The Copper grade distribution and spatial variability have been estimated through histograms (Figure 33) and variograms (Figure 35) computed on the well log data. We also denoted the correlation between Ag and Cu concentration (Figure 34) with a correlation coefficient of 0.67, although it has not been used in the simulation process because of the limited amount of data in the sandstone layer.

<figure><img src="image_90.png" /><figcaption>**Figure 33: Copper grade distribution in the Kupferschiefer sandstone layer based on WW well logs (log scale)**</figcaption></figure>

<figure><img src="image_91.png" /><figcaption>**Figure 34: Crossplot of Cu versus Ag in the Kupferschiefer deposit (in log-log scale)**</figcaption></figure>

These parameters have been used in a Sequential Gaussian Simulation process with an Ordinary Kriging that populated the Kupferschiefer Sandstone layer of the grid with Copper grades (Figure 36).

<figure><img src="image_92.png" /></figure>
<figure><img src="image_93.png" /><figcaption>Figure 35: Variogram for Copper in the Kupferschiefer sandstone</figcaption></figure>

<figure><img src="image_94.png" /></figure>
<figure><img src="image_95.png" /><figcaption>Figure 36: Simulated Copper grades in the Kupferschiefer sandstone in mass percent (log-scale)</figcaption></figure>
# 5. Hydrogeological model

## 5.1. Characterisation and simplification of aquifers

A hydrogeological model has been built on the basis of the simplified regional structural model (Section 2.5). This model aims at providing regional boundary conditions for further simulations (e.g. permeability stimulation as carried in WP2.3), and hydrogeological insights for environmental studies (as part of WP5). In this model, the fault network has been simplified to adapt its resolution to the scale of the study. The modelling choices that have led to the simplification are meant to limit their impact on the fluid flow behaviour. Major faults and principal directions have been kept.

The location of aquifers has been identified from the literature and from the study of the four available logs from Weisswasser region (Figure 18). Values of salinity, porosity and permeability have been estimated to evaluate the quality of each aquifer. The main aquifers are located from Buntsandstein to Quaternary formations (Figure 37).

It appears on Figure 18 that there are some aquifers in the Cainozoic formations (Tertiary and Quaternary) even if there is a lack of data for the Quaternary. Some aquifers, which are exploited for the consumption of the population, much have a low salinity. From the literature, Buntsandstein formations are known to have a high salinity, which is consistent with the logs on Figure 18 where the Buntsandstein formation is shown to be between two anhydrite layers. With this information we can infer that the drilling mud had a high salinity. Salinity in the wells has been evaluated relatively to the drilling mud. It appears that there is a decrease of salinity from Muschelkalk formation to Quaternary. So aquifers used for consumption is probably aquifers of Tertiary and Quaternary formations.

According to the study of logs and data of the bibliographic study, the aquifers network is reduced to two aquifers located from the Buntsandstein to Quaternary formations (Figure 37). At some places, these aquifers are connected because of the erosion of the corresponding formations which are bad aquifers (like Muschelkalk or Keuper).

## 5.2. Boundary conditions

Some boundary conditions such as the rate of rainfall and the rate of infiltration in the first meters are required for constraining the simulations. These data are available in the literature. The study area appears to be on a same watershed of the *Oder*.

According to the geological map, we know that aquifers reach the surface and thus are not confined. It means there is a possible direct water supply from the surface as shown in Figure 38. Based on the average precipitation rate in the area we would consider a water influx of about 500mm/year for the study area. To account for the actual water influx we would have to consider other factors which have an influence on the infiltration rate such as the type of soil, the kind of material that covers the ground, the soil moisture, the topography.

<figure><img src="image_96.png" /></figure>
<figure><img src="image_97.png" /><figcaption>**Figure 37: Identification of the various possible aquifers**</figcaption></figure>

<figure><img src="image_98.png" /><figcaption>**Figure 38: Cross section of aquifers and possible water flow**</figcaption></figure>
### 5.3. Preparations for fluid flows modelling

For the fluid flows modelling we have to build a grid in which we attribute some useful properties: Net to Gross, permeability, porosity, pore volume, cell volume. The construction of the grid is done in the "Flow Simulation Grid" workflow on SKUA. Then the attribution of properties is performed with "Reservoir Properties" workflow on SKUA.

<table><caption>Table 7: Data on aquifers of the study picked up in Poland from Tarkowski et al. (2009)</caption><thead><tr><td>Formation</td><td>Area (km²)</td><td>Porosity (%)</td><td>Net gross ratio (%)</td><td>Storage Capacity (Mt)</td></tr></thead><tbody><tr><td>Lower</td><td>24 562</td><td>20.5</td><td>40</td><td>7 646.9</td></tr><tr><td>Lower Trias</td><td>112 036</td><td>9.7</td><td>70</td><td>26 494.1</td></tr></tbody></table>

Because of the lack of data, the Net To Gross ratio (NTG) has been inferred for each stratigraphic unit thanks to the data found in bibliography (cf. Table 7). Here any cell with a connected porosity above 2% is considered as productive. The NTG for each stratigraphic unit has been estimated from this threshold (Table 8) knowing that it equals to the ratio between the volume of productive cells and the total volume of the stratigraphic layer.

<table><caption>Table 8: Estimation of NTG for each stratigraphic unit</caption><thead><tr><td>Stratigraphic unit</td><td>NTG</td></tr></thead><tbody><tr><td>Quaternary</td><td>0.88</td></tr><tr><td>Tertiary</td><td>0.4</td></tr><tr><td>Upper Cretaceous</td><td>0.74</td></tr><tr><td>Lower Cretaceous</td><td>0.4</td></tr><tr><td>Keuper</td><td>0.54</td></tr><tr><td>Muschelkalk</td><td>0.38</td></tr><tr><td>Buntsandstein</td><td>0.75</td></tr><tr><td>Aller</td><td>0.68</td></tr><tr><td>Leine</td><td>0.42</td></tr><tr><td>Strassfurt</td><td>0.24</td></tr><tr><td>Oberer Werra Anhydrite</td><td>0.24</td></tr><tr><td>Aeltestes Steinsalz</td><td>0.24</td></tr><tr><td>Unterer Werra Anhydrite</td><td>0.05</td></tr><tr><td>Zechsteindolomit to Rotliegendes</td><td>0</td></tr></tbody></table>

$$ \text{pore volume} = \text{NTG} * \text{cell volume} * \text{porosity} $$

The pore volume is computed as:

<figure><img src="image_99.png" /></figure>
The other necessary properties have been calculated thanks to SKUA-GOCAD tools.

5.4. Export of flow grid on Eclipse for the fluid flows modelling

The resulting grid and properties have been exported to Eclipse format (Figure 39).
We choose to export the objects in ASCII format to make it possible to manually edit
the files. ASCII export also provides more options than the binary export.

The grid is required to be built with straight pillars and it is easier to use stair steps
along unconformities. After building the grid with the GOCAD-SKUA workflow, the non-
neighbour connections are estimated to relate cells separated by faults or
unconformities. Inactive cells, aka dead cells, are encoded in the ACTNUM property,
in order to inactive dead cells in Eclipse for the simulation.

<figure><img src="image_100.png" /><figcaption>Figure 39: District model flow grid loaded in Eclipse flowviz.</figcaption></figure>

<figure><img src="image_101.png" /></figure>
6. Partner communication and model export

The stratigraphic, fault and boundary surfaces of the different structural models have
been exported as different standard or open files formats:

• .ts (Gocad TSurf format)

• .obj (http://en.wikipedia.org/wiki/Wavefront_obj_file)

* .mesh, .meshb (ASCII and binary, http://www-roc.inria.fr/gamma/gamma/Membres/CIPD/Loic.Marechal/Research/LM5.html)

* .ply (ASCII and binary, single and double precision, http://en.wikipedia.org/wiki/PLY_(file_format) )

• .off (http://www geomview org docs html OFF.html)

• .stl (ASCII and binary, http://en.wikipedia.org/wiki/STL_(file_format))

For each format, a folder containing a file per surface has been generated. The name
of the file reflects the type of surface it contains. The first part skua_model_
corresponds to the method for generating the original surface. Then, there is the type
of surface: fault_ for faults, horizon_ for geological horizons including stratigraphic
horizons and unconformities, or boundary_ for each surface of the bounding box of the
study volume. Finally, the last part corresponds to the name of the corresponding
surface in the stratigraphic column. The topographic surface is named
skua_model horizon Surface_ts.

The different surfaces have been exported from the Gocad-Skua model in .ts format.
A program for converting surface files has been created using RINGMesh library
(http://www.ring-team.org/index.php/software/ringmesh) and the Geogram tools it
embeds. The result can be visualised with MeshLab<fnref n="3" /> (Figure 40).

<figure><img src="image_102.png" /><figcaption>Figure 40: The final model in .obj format loaded in Meshlab software.</figcaption></figure>

<figure><img src="image_103.png" /></figure>

<footnote n="3">³ http://meshlab.sourceforge.net/</footnote>
## 7. Discussion and conclusions

Geological models for various purposes have been produced at different scales from the available data in the BIOMOre project:

* Detailed regional scale:

* Structural model (based on ProMine model)

* Simplified regional scale:

* Structural model

* Stratigraphic grid with petrophysical properties

* Hydrogeological model for Eclipse

* District scale:

* Structural model

* Stratigraphic grid with petrophysical properties

* Well field scale (preliminary):

* Regular grid

The different models have been exported in appropriate formats to make them easily usable by BIOMOre partners, i.e. without requiring SKUA-GOCAD when possible. Structural models are shared in open standard formats as a set of geological surfaces. This required developing a homemade executable based on RINGMesh library. Hydrogeological models are in Eclipse format as it is relatively standard for this kind of application and readily available as a Gocad export.

This suit of models gives a basis for the modelling strategy to be applied for further geological modelling tasks in the BIOMOre project and any following study for planning pilot plants or actual industrial scale plant.

<figure><img src="image_104.png" /></figure>
<figure><img src="image_105.png" /><figcaption>Appendix 1. Cross sections</figcaption></figure>
<figure><img src="image_106.png" /><figcaption>Figure 41: Cross section number 1 without vertical exaggeration</figcaption></figure>
## Appendix 2. Well field scale models

A short study has been carried to present quick rock stimulation indicators, which could be used for preparing well-design. It is presented in this appendix as an example of complementary tool for rapid behaviour appraisal. It should be noted that this approach may only be used as a preliminary evaluation and that more evolved rock stimulation simulations (such as carried in the task 2.3) are required for a thorough study.

### 1.1. Modelling the well field design for ISL

The ISL wells design has to take into account several aspects, including geology, hydrogeology, rock mechanics and operational aspects. However, this study does not aim to make technical specification of drilling methods to be used, nor a comprehensive economical appraisal, but it is an attempt to compare two designs from technical and economical point of view. For this, one can calculate capital expenditures to mount an installation required for one well and possible revenue from copper extraction, depending of the leaching area. In this case, the characteristics of the well location (WW1103) and two well patterns (radial and rectangular) are used. The parameters of wells (diameter, length) are taken from the similar projects for shale gas production (Table 9).

Combination of vertical and horizontal wells are planned to have greater accessibility within tabular ore body. The wells assumed to be hydraulically fractured in several stages. The technique to be used in the HF assumes that two wells are alternately fracked ("zipperfrac"). It should be noted that the average thickness of shale strata at Weisswasser is less than 1 m that can make vertical drilling complicated. The patterns will be used for simulating the fracture propagation using the Gocad/SKUA Wizard plugin.

#### 1.1.1. Radial pattern

The radial pattern assumes that twelve (12) wells are drilled from the pad down to the Kupferschiefer strata, then inclined to become horizontal and to have symmetrical radial dispersion within the shale formation. The lateral horizontal extend of the wells is about 800 m, creating a circular drained area. Injecting and production wells are assumed to alternate laterally. A typical well field plan using horizontal wells is shown on Figure 42.

<figure><img src="image_107.png" /></figure>
<figure><img src="image_108.png" /><figcaption>**Figure 42: Radial well field pattern**
comprising 12 horizontal wells drilled within the Kupferschiefer shale formation.</figcaption></figure>

The wells were built in Gocad/SKUA by drawing vertical and horizontal trajectory curves from which the wells were created. The horizontal curves were drawn on the Kupferschiefer horizon to follow strictly the shale formation with a lateral extend of 800m and were move upward by 0.5 m to be within the formation. In vicinity of the NNE-SSW trending fault, they were stopped due to the down-throw along the fault.

#### 1.1.2. Rectangular pattern

In the rectangular pattern case, eleven (11) horizontal wells are assumed to be drilled from two vertical wells (pads) – 5 injection wells and 6 production wells. The wells split in such a way to create ‘zip’ structure (injection wells are surrounded by the production) within the Kupferschiefer shale, forming a rectangular area in plan. Since the study area is slightly dipping, the injection wells are located at slightly higher altitude than production wells to take advantage of the gravity. Each vertical well is 800 m long and separated by 150 m span from the other wells (Vermylen, 2011).

<figure><img src="image_109.png" /><figcaption>**Figure 43: Rectangular well field pattern**
comprising five injection and six production wells drilled from two pads within the Kupferschiefer shale formation.</figcaption></figure>

The geometrical parameters are used for calculating the costs and revenues in the preliminary economic appraisal. The total costs can be estimated by the length of

<figure><img src="image_110.png" /></figure>
vertical and horizontal drilling for one well and amount of copper extracted from the covered area. The input data for calculation are summarized in Figure 43.

<table><caption>Table 9: Summary of characteristics of both radial and vertical well field design.</caption><thead><tr><td></td><td>Number of wells</td><td>Horizontal length, m</td><td>Vertical length, m</td><td>Covered area, km²</td><td>Spacing between wells</td></tr></thead><tbody><tr><td>Radial</td><td>12</td><td>1,725</td><td>800</td><td>2.0</td><td>30</td></tr><tr><td>Rectangular</td><td>11</td><td>3,460</td><td>800</td><td>1.2</td><td>150</td></tr></tbody></table>

<figure><img src="image_111.png" /></figure>
## 1.2. Geomechanical characterisation

Intrinsic parameters (strength of the rock, angle of internal friction) and in-situ stress field control the propagation of hydraulic fractures. Some of geomechanical parameters required for HF simulation were taken from the KGHM database for shales (Table 10); additional unknown parameters are taken from the published literature.

<table><caption>Table 10: Geomechanical properties of shale used in hydraulic fracturing simulation.</caption><thead><tr><td>Rock tensile strength T (MPa)</td><td>Poisson's ratio v</td><td>Internal friction angle μ</td><td>Young's Modulus E (GPa)</td></tr></thead><tbody><tr><td>5</td><td>0.4</td><td>0.36</td><td>35</td></tr></tbody></table>

The measurements of stress magnitudes in the Weisswasser region have not been conducted (or are not available to our knowledge). Ongoing investigation of the stress regime in the study area conducted by (Konietzky & Weber, 2016) revealed four stress measurements around the study area. The data shows presence of all three fault regime: two points with strike-slip (SS) regime, one in compressive (reverse RS) and one in extensional (normal NR), following the Anderson's classification of the stress (Zoback 2007). These controversial observations suggest that there have been change in the stress field over the time.

There is large correspondence between the in-situ stress field and large-scale active faults in the area (Zoback 2007). Looking on the character of the faults in the study area, it deems that the Weisswasser region has a normal tectonic regime, as it contains mainly normal faults. In a normal faulting environment, it is assumed that $S_v > SH_{max} > Sh_{min}$. Vertical stress magnitude can be found as the weight of the overburden (Zoback, 2007):

$$S_v = \int_{0}^{z} \rho(z) g dz;$$

Where $S_v$ is vertical stress, $\rho(z)$ is the density at depth $z$ and $g$ is gravity (9.81 m.s$^{-2}$). Therefore, the vertical stress increases with depth at a rate of about 23 MPa.km$^{-1}$ (Zoback, 2007). The magnitudes of the two horizontal stresses can be found by the expression:

$$S_H = \frac{\nu}{1 - \nu} S_v;$$

where $\nu$ is the Poisson's ratio. In the normal faulting regime, if the pore pressure is close to hydrostatic one, the least principal stress can be estimated as 0.6 $S_v$ (Zoback, 2007) assuming a Poisson coefficient at 0.4.

Pore pressure is created by the fluid phase trapped in the pore space of the rock. Usually it can be estimated as the drilling mud weight (Zoback, 2007); however, no data is available. For this particular case the pore pressure is assumed to be hydrostatic with a

<figure><img src="image_112.png" /></figure>
rate of 10 MPa.km⁻¹. In fact, the pore pressure can be higher, because the shale formation located below the anhydrite layer that could cause overpressure conditions by trapping and confining the fluids.

### 1.3. HF simulation

The HF Wizard Gocad plugin was developed by Chaumont & Cosson, (2012) for prediction of fracture occurrence and orientation against the distance from a injection well. The Wizard is created for both isotopic and anisotropic, but in this work, only the isotropic case was studied. The anisotropic case requests the evaluation of the strain tensor which is unknown.

The Wizard works with a Voxet in which the required parameters are stored in each cell. The HF simulation is modelled only for the Kupferschiefer shales, i.e. it contains exclusively the properties of shales. For each well field design the Voxets were created from the object cage with two regions representing the wells and the surrounding rock. The sizes of the Voxets for radial pattern are 2140 m × 1850 m × 575 m and for rectangular is 2040 m × 1080 m × 350 m. It should be mentioned that the Voxet dimension is far from the shale layer in geological grid, but it was required for the modelling purposes. The Wizard allows the evaluation of the required rock properties at each cell node with a non-conditional Gaussian simulation (SGS). Thus, the following parameters are required:

* Rock tensile strength (T_simulation);

* Poisson's coefficient (poisson_simulation);

* Internal friction angle φ (phy);

* Rigidity μ (mu).

The hydraulic pressure created by fluid injection dissipates with an increasing distance r from the well (Figure 44). To model the change in injection pressure at a radial distance r from the well, the Dupuit-Forchheimer formula is used (Royer, 2014):

$$P(r) = -P_{inj} \left( \frac{\ln(r_0)}{\ln(r)} \right);$$

where P_inj is the injected water pressure (MPa), r₀ is the diameter of the well and r is radial distance from the well. The radial distance is computed using the 'neighbor distance' Gocad function from the well region. Note here that for the rectangular pattern, the radial distance is uniform along the well, while for the radial pattern the difference in radial distance is increasing that can create stand-alone areas on the periphery. It can be dangerous for the ISL performance since this areas are too distant from well and can become unfractured.

<figure><img src="image_113.png" /></figure>
<figure><img src="image_114.png" /><figcaption>Figure 44: Injection pressure propagation in the rectangular (left) and radial (right) wellfield designs.</figcaption></figure>

Then, all the rock properties are simulated at each cell node in the Voxet, the other constant parameters are required for computation (Table 11). Regarding the choice of the parameters, the volumetric mass is an average density of the rock (not shale density) since it is used for determining the horizontal stress; the reference altitude is the z coordinate of surface; $σ_1$ is the tectonic stress which is equal to the difference between the minimum $S_h$ and maximum $S_H$ horizontal principal stresses. For the isotopic case, it is assumed that $S_h = S_H$, hence the tectonic stress equals to zero. The results of the parameter simulation are summarized in Table 11.

<table><caption>Table 11: Parameters used for hydraulic fracturing simulation.</caption><thead><tr><td>Parameters</td><td>Units</td><td>Values</td></tr></thead><tbody><tr><td>Volumetric mass</td><td>kg.m<sup>-3</sup></td><td>2,650</td></tr><tr><td>Well diameter,</td><td>m</td><td>0.21</td></tr><tr><td>Rock tensile strength</td><td>MPa</td><td>5</td></tr><tr><td>Water initial pressure</td><td>Mpa</td><td>5; 10; 15</td></tr><tr><td>Reference altitude,</td><td>m</td><td>100</td></tr><tr><td>Reference pressure</td><td>MPa</td><td>0.017</td></tr><tr><td>sigma_1, sigma_2, sigma_3</td><td>MPa</td><td>0</td></tr></tbody></table>

It is important to describe the basic outline of calculation in the Wizard. The script calculates firstly the principal stresses and normalizes them to the tensile strength. Next, it determines the tectonic regime (normal, strike-slip or reverse) and calculates the geometrical parameters of the Mohr circles (center, radius). Therefore, this is the parameters determined by the stress field that acts in the rock. The critical domain can be found through the intrinsic parameters of the rock (angle of internal friction, cohesion and tensile strength), depending on the criterion applied (straight line for the Coulomb domain and parabola for the Griffith). The fracturing will occur in case the water injection pressure $P_f$ exceeds the resulting normal pressure $σ_{ex}$, where the resulting normal

<figure><img src="image_115.png" /></figure>
pressure is the difference between the center of the Mohr circle and the normal stress at failure plane (Figure 45).

<figure><img src="image_116.png" /><figcaption>**Figure 45: Visualisation of failure criterion (after Chaumont & Cosson, 2012).**</figcaption></figure>

**1.4. Results and discussions**

As expected, the fracture propagation depends on the initial water pressure for every pattern. Figure 46 shows the fracture evolution depending on the magnitude of initial water pressure for the rectangular pattern. The distribution of fractures is symmetrical in a way that allows the determination of the critical pressure at which the fracture fronts get connected between the wells. For the rectangular pattern, the critical pressure is around 13.5 MPa. Hence, the advantage of this pattern is in relative predictability of fracture propagation.

Fracture propagation in the radial pattern can be seen in Figure 47 for injecting water pressures of 5, 10 and 15 MPa. It is noticeable that the fracture propagation is less uniform and at $P_{inj}$ = 15 MPa the most remote areas is not as intensively fracked as for the central part. It can be explained by the varying difference in the radial distance between wells. For the same values of injection pressure, the resulting fracture network is less connected compared with the rectangular pattern. The critical pressure at which the fracture front get connected between wells depends on the length of vertical well and the angle between wells (or number of wells).

Consequently, it can create a severe impact on the lixiviant flow during leaching phase: the production on remote areas can be hindered since some patches of the shales on the periphery will not be subjected to the fracking. Therefore, it is more complicated to predict the behaviour of hydraulic fracturing for the radial pattern.

<figure><img src="image_117.png" /></figure>
<figure><img src="image_118.png" /></figure>

<figure><img src="image_119.png" /></figure>

<figure><img src="image_120.png" /></figure>

**Figure 46: Hydraulic fracture propagation**
(in red) for the rectangular wellfield for injection pressure of 5 MPa (A), 10 MPa (B) and 15 MPa (C).

<figure><img src="image_121.png" /></figure>
A

<figure><img src="image_122.png" /></figure>

B

<figure><img src="image_123.png" /></figure>

C

<figure><img src="image_124.png" /></figure>

Figure 47: Hydraulic fracture propagation in the radial wellfield for injection pressure of 5 MPa (A), 10 MPa (B) and 15 MPa (C).

<figure><img src="image_125.png" /></figure>
1.5. Conclusions

In order to provide a preliminary investigation of fracture propagation in the shale formation, two well patterns have been studied: rectangular and radial. The synthetic block models were created in the Gocad/SKUA for performing the HF simulation, developed by Chaumont & Cosson, (2012). The parameters for calculation were updated to the present case.

The model is based on the assumption that the study area has normal faulting regime. The magnitudes of principal stresses are estimated as S<sub>v</sub> = 23 MPa.km<sup>-1</sup>, S<sub>H</sub> = S<sub>h</sub> = 13.8 MPa.km<sup>-1</sup> and fracture propagation is defined by the Coulomb-Mohr criterion. The intrinsic parameters were taken for the samples from the Rudna mine.

It was found that the rectangular pattern produces more uniform fracturing along the wells, while there is considerable variability in generating fractures in the radial pattern for the same length of wells (800 m). The critical pressure required to create interconnected fracture network is about 13.5 MPa for the rectangular pattern. However, it is not trivial to estimate it for the radial pattern. Indeed, the radial pattern may reveal the presence of 'blind areas' where the fractures do not fully connect between wells and the excessive fracturing in the vicinity of the central well, which can propagate to the adjacent layers (dolomites and sandstones). On the other hand, the blind areas can serve as pillars in the same manner as it in room and pillar mining method. Thus, in case of dissolution of carbonate-cemented rock fabric mineral, there will be no collapse of fracked volumes and as a consequence, no surface subsidence. However, the presence of the pillars automatically means loss in metal extraction.

Additionally, the direction of fractures created and/or reactivated dependent on interconnection of physical properties of shales, in-situ tectonic stress, number of existing fractures and pore pressure. It has to be investigated to avoid acid leakage and ground water contamination.

To conclude, the key geological/geomechanical issues that affect the modelling of hydraulic fracturing can be summarized as follows:

* Limited understanding of the *in-situ stress* in the study area is resulting in inability to correctly predict fracture patterns that have importance for environmental study;

• Identification of the shale formation boundaries and geometry;

* Lack of data about pre-existing fracture network.

<figure><img src="image_126.png" /></figure>