Digital rock physics approach to simulate hydraulic effects of anhydrite cement in Bentheim sandstone

Cementation of potential reservoir rocks is a geological risk, which may strongly reduce the productivity and injectivity of a reservoir, and hence prevent utilisation of the geologic subsurface, as it was the case for the geothermal well of Allermöhe, Germany. Several field, laboratory and numerical studies examined the observed anhydrite cementation to understand the underlying processes and permeability evolution of the sandstone. In the present study, a digital rock physics approach is used to calculate the permeability variation of a highly resolved three-dimensional model of a Bentheim sandstone. Porosity-permeability relations are determined for reactionand transport-controlled precipitation regimes, whereby the experimentally observed strong decrease in permeability can be approximated by the transport-limited precipitation assuming mineral growth in regions of high flow velocities. It is characterised by a predominant clogging of pore throats, resulting in a drastic reduction in connectivity of the pore network and can be quantified by a power law with an exponent above ten. Since the location of precipitation within the pore space is crucial for the hydraulic rock properties at the macro scale, the determined porosity-permeability relations should be accounted for in large-scale numerical simulation models to improve their predictive capabilities.


Introduction
Mineral dissolution and precipitation are micro-scale processes, which may significantly alter the internal rock structure and consequently affect the effective hydraulic behaviour of the system at the macro-scale. Predicting changes in rock permeability is of paramount importance for most applications related to geologic subsurface utilisation (Jacquey et al., 2015;Kleinitz et al., 2001;Regenspurg et al., 2015), especially regarding the productivity and injectivity of a reservoir. Permeability is without question the most crucial hydrogeologic parameter, but it is often difficult to determine due to its enormous variation over space and time in natural geological systems (Hommel et al., 2018;Bernabé et al., 2003). Usually, there is a strong positive correlation between porosity and permeability, which can be simply quantified by an exponential equation related to the initial values of porosity and permeability. The relevant exponent varies depending on the represented particular geochemical reaction or hydrothermal reservoir system (Kühn, 2009).
For geothermal exploration, a drill hole was deepened at Allermöhe (Hamburg, Germany) into the Rhaetian sandstones. The originally open pore spaces were found to be filled with anhydrite to a large extent (Baermann et al., 2000a). Therefore, the potential for geothermal exploitation was limited and the drill hole had to be abandoned. Beforehand, extensive studies were performed to reveal insights into the regional distribution of the cementation (Baermann et al., 2000b). Pape et al. (2005) analysed the observed local clogging phenomena due to anhydrite precipitation in the porespace of related drill hole samples and studied the pore size geometry. They distinguished between two facies types with different diagenetic history: (1) a fine-grained and uncemented sand with small pores, mechanically compacted during diagenesis and (2) a coarser sand with larger pores, almost completely cemented by anhydrite after its compaction.
In order to understand the observed anhydrite cementation in the core samples from the borehole at Allermöhe,  performed fully coupled numerical simulations of reactive flow at the local and regional scales. A reser- Figure 1. Binarised micro-CT scan of the Bentheim sandstone has an initial porosity of 23 % and comprises pore space (blue) and grains (brown). Exemplary slices through the digital rock show the anhydrite cementation (yellow) for a reaction-and a transport-controlled precipitation at a porosity of 17 %.
voir model was used to simulate the chemical stimulation of the sandstone formation via a forced increase in porosity and permeability in the close borehole vicinity. At the reservoir scale, numerical simulations of Kühn and Günther (2007) indicate that strata-bound convective flow in the Rhaetian sandstone reservoir is insufficient to explain the observed high degree of anhydrite cementation. However, it is shown by  that brine is forced to precipitate around fracture zones, if hot fluids flow up along faults, heating up the aquifer.
A core flooding experiment simulating the conditions in an operated geothermal reservoir was performed to validate the abovementioned numerical models, applied for the Allermöhe site (Bartels et al., 2002). A Bentheim sandstone sample was prepared with anhydrite and subjected to a temperature gradient, while the changes in permeability were measured along the core. Within the cold upstream and hot downstream regions of the core, anhydrite was dissolved and precipitated, respectively.
Within the presented study, we apply a digital rock physics approach to determine if observed anhydrite precipitation in the core flooding experiment of Bartels et al. (2002) and samples of the Allermöhe site (Baermann et al., 2000a;Pape et al., 2005) can be ascribed in particular to a reactionor transport-controlled process. To upscale the chemical reactions from the micro to the macro scale, we determine the required porosity-permeability relationship. The impact of precipitation on the effective hydraulic rock properties is especially of interest for a sustainable utilisation of the geological subsurface and for risk assessment to avoid unsuccessful drilling campaigns like the one at Allermöhe.

Material and methods
A digital sample of the Bentheim sandstone is used to determine the dominant mechanism of pore space alteration due to anhydrite cementation, since its mechanical and hydraulic properties have been well investigated by several studies (Bakker and Barnhoorn, 2019;Jacquey et al., 2016;Peksa et al., 2015;Wetzel et al., 2018) and it represents a common reference reservoir rock. The microstructure of the Bentheim sandstone derives from a subset of a publicly available micro-CT scan (Herring et al., 2018 and comprises 500 3 voxels with a resolution of 4.9 µm (Fig. 1). The permeability of the binarised 3D image is calculated from the flow field by solving the steady-state Stokes equation by means of the OpenFOAM software package (Weller et al., 1998). The calculation of the pore network is done in Python using the PoreSpy package (Gostick et al., 2019). The investigated digital sample has a porosity of 23.4 %, which is comparable to the sandstone of the core flooding experiments used by Bartels et al. (2002) with 24.6 %. They simulated conditions in a managed geothermal reservoir and investigated the effect of anhydrite dissolution and precipitation on permeability.
The present digital rock physics approach assumes a flow velocity dependency of mineral precipitation to estimate the porosity-permeability relation instead of performing fullycoupled reactive transport simulations. The voxel-based algorithm converts pore space to minerals, depending on the magnitude of the flow velocity, since it reflects the local pore morphology as well as the connectivity of the entire pore space. Two contrasting scenarios are simulated associated with the general geochemical reaction regimes: (1) a preferential precipitation at high flow velocities (> 75th percentile), which can be associated with a transport-controlled process where the reaction is limited by the availability of re-  actants transported to the fluid-mineral interface (Noiriel and Daval, 2017), and (2) a uniform precipitation independent of the fluid flow, that represents a surface reaction-controlled process, whereby the chemical reaction rate limits precipitation and pore space alters homogeneously (Fig. 1). A detailed description of the model set-up as well as the iterative precipitation approach is described in Wetzel et al. (2020).

Results
The geochemical reaction regime strongly impacts the evolution in permeability of the digital rock sample as illustrated by the large variations between transport-controlled and reaction-controlled porosity-permeability relations. In case of precipitation at high flow velocities, a permeability decrease of one order of magnitude occurs for a reduction in porosity to 20 %, whereas in the case of a uniform pore alteration this permeability is achieved at a substantially lower porosity of 13.5 % (Fig. 2a). The flow velocity thresholds represent the intensity of the dominant mechanism regarding both contrasting processes and allow for gradation between them. Permeability evolution due to anhydrite cementation in a Bentheim sandstone with comparable initial hydraulic properties as for the chosen digital sample is described by Bartels et al. (2002, Table 1). The observed strong decrease in permeability can be approximated by the simulated porositypermeability relation for a transport-controlled precipitation, assuming mineral growth in regions comprising flow velocities > 75th percentile (Fig. 2b). Further, the strong decrease in permeability can be described by a power law with a high exponent (η) above 11.3 (Bartels et al., 2002). This is similar to the findings of Baermann et al. (2000a, b), who determined an exponent of around twelve from field samples taken at the Allermöhe well.
Since the experimentally observed permeability changes due to anhydrite cementation can be associated with the simulated transport-controlled precipitation, the digital rock physics approach further allows for quantification of significant alterations in the pore morphology. Noticeable changes in the number and sizes of throats indicate a preferential clogging of pore throats. Compared to the unaffected microstructure, 23 % of the initial pore throats are completely closed at a reduced porosity of 17 %. Moreover, the median throat size substantially decreases from 21 to 11 µm, while particularly the diameter of larger throats is also reduced (Fig. 3a). The mean pore diameter reduces only slightly from 92 to 82 µm. The pore size distribution becomes more right-skewed and the number of smaller pores considerably increases (Fig. 3b), whereas the absolute number of pores does not change due to precipitation. The reduction in pore diameters at a reduced porosity of 17 % exhibits a broad right-skewed distribution with a median reduction in the pore diameter of 10 µm. The diagram presented in Fig. 3c, which relates the changes in pore sizes to the initial diameter, demonstrates comparably lower reductions in the diameters of smaller pores than in those of the larger pores, while none of the pores are completely closed.
The initial flow regime within the pore space of the Bentheim sandstone is highly variable. There are well interconnected regions with preferential flow paths and comparably high flow velocities as illustrated by the streamlines in Fig. 4. The simulated transport-controlled pore alteration strongly impacts the spatial distribution of fluid flow paths and velocities. At a reduced porosity of 21 %, flow velocities noticeably decrease. Successive precipitation leads to the progressive closure of flow paths due to the described preferential precipitation in the pore throats. At a porosity of 17 %, the majority of initial flow paths are cut-off, and the pore space is considerably less connected (Fig. 4). These results are in accordance with the different facies types found by Pape et al. (2005), but also support the conclusion of Bartels et al. (2000) on the permeability of the Rhaetian sandstone at Allermöhe being particularly reduced due to the narrowing of pore throats as a result of diagenetic anhydrite precipitation.

Discussion and conclusions
Mining for geothermal heat usually requires at least one or more production and injection wells. These installations are the major capital investments and demand economic feasibility, which mainly depends on the permeability and porosity of the reservoir. A major geological risk is the clogging of the pore space of the reservoir rocks by precipitation of minerals such as anhydrite , as it was the case for the deepened Allermöhe, where the Rhaetian sandstone was largely cemented with anhydrite. Consequently, determination of the origin of the anhydrite cement pore fillings represents a basis for exploration to lower the risk for failure for future projects. This is mainly investigated by fully coupled numerical simulations of flow, heat transfer, species transport and chemical reactions at the reservoir scale (Kühn and Günther, 2007;. Precipitation reduces porosity in general, but its effect on transport properties strongly depends on its location. Various studies discuss if precipitation preferably occurs in larger pores or in contrast in smaller ones (Stack, 2015;Steinwinder and Beckingham, 2019). For the Allermöhe sandstone, it is presumed that larger pores were favoured for anhydrite precipitation during its diagenesis . This was also studied by means of a one-dimensional numerical test case carried out by Mürmann et al. (2013), who explained the heterogeneous anhydrite cementation patterns by a variation of the fluid supersaturation and fluid velocity within the pores.
In this study, we present similar results by applying a digital rock physics approach (Wetzel et al., 2020) with the advantage of a significantly higher resolution and representation of the specific porous network of the Bentheim sandstone. Thus, it is possible to adequately predict hydraulic rock properties by characterising the micro structure, morphology and connectivity of pores. In view of a virtual laboratory, digital rock models allows to consider varying distributions of secondary minerals for the same rock sample and can be used to establish trends of evolving geophysical properties, e.g., as consequence of diagenetic processes (Hosa and Wood, 2020;Singh et al., 2020). The results of this study indicate that particularly larger pores are reduced in size by mineral precipitation, as demonstrated by Mürmann et al. (2013). Even though all pore sizes are affected by precipitation, we observe that the amount of smaller pores considerably increases. Hence, the reduction in pore size is limited to a certain size, since pores are successively detached from the high velocity flow field and hence an assured reactant replenishment. Much more important and cru- Figure 4. Evolution of the streamlines for flow in x direction (red: high velocities, blue: low velocities) and the pore network (red: high pore diameter, blue: low pore diameter) considering a transport-controlled reaction regime with preferential precipitation in regions of high flow velocities (> 75th percentile). cial for the permeability development of the Bentheim sandstone are the pore throats. Their clogging is responsible for the steep porosity-permeability relation with an exponent of well above ten.
Of course, mineral precipitation is complex and not exclusively depending on flow velocity. The chemical reaction regime is controlled by various other factors like fluid chemistry, temperature, pressure, mineralogy, pore morphology and transport properties, whereby the impact of each aspect depends on the particular process. Moreover, the presented method bases on a micro-scale representation of the reservoir rock and therefore does not take larger scale heterogeneities into account (e.g. fissures, fractures or facies changes). Nevertheless, it proves to be a simple way to estimate porositypermeability relations conditional to the geochemical reaction regime without the need to apply complex and computationally expensive reactive transport simulations.
We conclude that the applied digital rock physics approach (Wetzel et al., 2020) is a viable model to describe the transport-controlled precipitation of anhydrite in Bentheim sandstones at the micro-scale and quantify changes of the effective porosity and permeability of the rock. Our results coincide with laboratory studies (Bartels et al., 2002), conducted to simulate field observations to reveal insights about the regional distribution of the cementation of the Rhaetian sandstone by anhydrite. The impact of the chemical precipitation regime on the effective hydraulic rock properties can be quantified by a porosity-permeability relation using a dig-ital rock physics approach. These findings can be now applied for the assessment of equivalent reservoirs to succeed with future projects aiming at utilisation of the geological subsurface in the North German Basin.
Author contributions. The three authors have equally contributed to this paper. MW and MK conceived and designed the simulations; MW performed the research; MW, TK and MK analysed the data; MW, MK and TK wrote the paper. All authors read and agreed to the published version of the manuscript.
Competing interests. The authors declare that they have no conflict of interest.
Special issue statement. This article is part of the special issue "European Geosciences Union General Assembly 2020, EGU Division Energy, Resources & Environment (ERE)". It is a result of the EGU General Assembly 2020, 4-8 May 2020.
Acknowledgements. We are grateful for the editorial handling and the constructive comments provided by the two anonymous reviewers.
Financial support. This research has been supported by the German Federal Ministry for Education and Research (BMBF) (Geo:N project GEOSMART (grant no. 03G0867E)).
The article processing charges for this open-access publication were covered by a Research Centre of the Helmholtz Association.
Review statement. This paper was edited by Johannes Miocic and reviewed by two anonymous referees.