Permeability and Mineralogy of the Újfalu Formation, Hungary, from Production Tests and Experimental Rock Characterization: Implications for Geothermal Heat Projects

Abstract

Hundreds of geothermal wells have been drilled in Hungary to exploit Pannonian Basin sandstones for district heating, agriculture, and industrial heating projects. Most of these sites suffer from reinjection issues, limiting efficient use of this vast geothermal resource and imposing significant extra costs for the required frequent workovers and maintenance. To better understand the cause of this issue requires details of reservoir rock porosity, permeability, and mineralogy. However, publicly available data for the properties of reservoir rocks at geothermal project sites in Hungary is typically very limited, because these projects often omit or limit data acquisition. Many hydrocarbon wells in the same rocks are more extensively documented, but their core, log, or production data are typically decades old and unavailable in the public domain. Furthermore, because many Pannonian sandstone formations are poorly consolidated, coring was always limited and the collected core often unsuitable for conventional analysis, only small remnant fragments typically being available from legacy hydrocarbon wells. This study aims to reduce this data gap and to showcase methods to derive reservoir properties without using core for flow experiments. The methods are thin-section analysis, XRD analysis and mercury intrusion porosimetry, and X-CT scanning followed by numerical flow simulation. We validate our results using permeability data from conventional production testing, demonstrating the effectiveness of our method for detailed reservoir characterization and to better constrain the lateral variation in reservoir properties across the Pannonian Basin. By eliminating the need for expensive bespoke coring to obtain reservoir properties, such analysis will contribute to reducing the capital cost of developing geothermal energy projects, thus facilitating decarbonization of global energy supply.

1. Introduction

Hundreds of geothermal wells have been drilled in the Pannonian Basin, Hungary, in recent decades for direct-use applications of geothermal heat, with more than 900 active wells producing ~10.7 PJ of heat in 2019, the main uses being for balneology (~34%), agriculture (~27%), and district heating (~21%) [1]. This extensive geothermal development reflects the high typical heat flow (~100 mW m⁻²) and geothermal gradient (~45 °C km⁻¹) [2] caused by the crustal thinning that accompanied lithospheric extension in the Early-Middle Miocene [3,4]. A key geothermal aquifer target in the Pannonian Basin is the sequence of deltaic and turbiditic sandstones deposited during the Pannonian stage of the Late Miocene (e.g., [5,6]). These rocks were emplaced at a time of dramatic environmental change, resulting in rapidly-changing sedimentary environments: the ancestral southward drainage of the Pannonian Basin to the Aegean Sea became disrupted, creating an endoreic lacustrine environment, before the modern drainage to the Black Sea via the River Danube became established [7,8,9]. Infill of the Pannonian Basin ranged from lacustrine turbidites of the Szolnok Fm. and the Algyö Fm., deltaic sandstones of the Újfalu Fm. (formerly known as the Törtel Fm. [10]), and fluvial sandstones of the Zagyva Fm., which form both prolific hydrocarbon reservoirs and also geothermal resources [10,11,12,13]. Strictly speaking, these named ‘formations’ are lithofacies associations with diachronous boundaries [14], hence the ‘interleaved’ boundaries in Figure 1.

Our study was initiated as part of an attempt to understand the low injectivity at the Mezőberény geothermal project site in SE Hungary, located in the Békés Basin, one of the deep sub-basins of the wider Pannonian Basin (Figure 1). Given the delta plain depositional environment of the Újfalu Fm., this formation consists of a mix of sandstone, representing deposition within former channels, and siltstone/mudstone, representing overbank deposition (e.g., [16,17]). The lack of continuity of these sand bodies in any particular direction (for example, in the direction between the two wells at Mezőberény in Figure 1A) potentially has a clear adverse effect on reservoir properties [18].

Despite the long experience of production of geothermal heat and hydrocarbons from these sandstones, many geothermal sites experience reinjection issues; as a result, less than 10% of all Hungarian geothermal wells have been used as reinjection wells in the past decades (e.g., [5,6]). Chemical factors, such as scaling and clay mobilization, have previously been discussed as causes of these injection problems (e.g., [19]), but effective mitigation has not been developed. These injection issues have significant adverse financial impact on projects by limiting production rates and causing high workover costs, even forcing closure of some projects. A major hurdle for solving these problems has been the limited availability of published information on these issues and the reservoir properties in general. Subsurface data from geothermal sites is limited because these are typically managed and designed by small companies and municipalities, which cannot afford investment in subsurface data collection, as is standard for oil and gas companies and crucial for optimization of well planning and operations. In May 2021 the Mining and Geological Survey of Hungary (Magyar Bányászati és Földtani SZolgálat; MBFSZ) launched their web portal providing access to geological data, including subsurface data [20,21]. However, this does not address the issue of site-specific data, relevant to the present study, being unavailable. One of few recently published studies, by Varga et al. [22], on reservoir properties of Pannonian sandstone, provides a qualitative evaluation of scaling and cementation for a geothermal site at Szeged (Figure 1). In contrast, one of the most recent publications to quantify porosity and permeability of Pannonian sandstones dates from 1994 [11]. Istan Almási’s Ph.D. thesis, examined in 2001 [23], which reports analyses of 12,899 core samples from 689 petroleum wells in Hungary, is widely cited as a source of information on hydraulic properties of geothermal reservoir rocks. However, this thesis did not present any original experimental data, it is essentially a compilation of pre-existing data from diverse sources (e.g., [24,25,26,27,28,29]), most of which are inaccessible. This thesis, which is available online [23], reports summaries of results, the data compilation being provided only on a CD-ROM, which was unavailable to us; this dataset, held by MOL, remains confidential. Furthermore, this thesis presents aggregated results for the various sandstone facies present, there being no differentiation between the Szolnok Fm. (turbiditic), Újfalu Fm. (deltaic), and Zagyva Fm. (fluvial); as Figure 1 indicates, core samples from circa 2 km depth might be from any of these subdivisions. Nonetheless, Almási [23] reported that for reservoir sandstones from ~2 km depth the horizontal permeability is typically ~50% larger than the vertical permeability, both values being typically in the range ~100–200 mD; the porosity of these rocks is typically ~0.22 [23].

A major reason for the lack of recent quantitative studies of reservoir properties is that much of the core data was recovered decades ago and the amount of high-quality material that remains in the public domain is relatively limited. Furthermore, coring had low recovery because of the unconsolidated nature of many Pannonian sandstone formations. At the start of this study in 2018, we consulted MBFSZ and MOL Plc. (Magyar OLaj- és Gázipari Részvénytársaság, Hungarian Oil and Gas Public Limited Company), who maintain in Szolnok a core repository for petroleum wells in the Pannonian Basin. Access to core in the Újfalu Fm. from wells in the Békés Basin (Figure 1) was requested. Many core samples were available from the Újfalu Fm. elsewhere in Hungary and/or from other formations in the Békés Basin. However, for the Újfalu Fm. in the Békés Basin, the few core samples of sandstone that were collected, and were large enough for conventional laboratory experiments (i.e., plugging and flow experiments) to measure porosity and permeability, appear to have been processed long ago, the left-over fragments from wells in this area that were available to us from the MOL core repository being too small to be analyzed in this way.

This paper aims to begin filling this data gap on reservoir properties. We explore methods to extract porosity and permeability data from core material that is unsuitable for laboratory flow experiments. We obtained core for the Újfalu Fm. from the legacy Kond-1 and Gyoma-1 petroleum wells in the Békés Basin (Figure 1). These fragments are analyzed using thin section analysis, Mercury Intrusion Porosimetry (MIP), X-ray Diffraction (XRD), X-ray Computed Tomography (X-CT) imaging, and numerical flow simulations. The thin section analysis and XRD provide information on mineralogy, for the assessment of scaling and reservoir erosion risks. MIP, X-CT imaging, and numerical flow experiments provide new porosity and permeability data. We then compare the permeability values from these experiments with estimates from production tests from two geothermal wells at Mezőberény, also in the Békés Basin (Figure 1). These wells were drilled in 2012 to depths of ~2000 m and target the Újfalu Fm. [15,18]. With this approach we not only provide new data but also highlight the need for more reservoir characterization of Pannonian Sandstones. This is crucial for future efforts to avoid and solve injection problems and optimize exploitation of the enormous geothermal resource in this region. This paper is structured as follows: in section two our dataset is presented and the production tests and associated permeability estimates are discussed. In section three we describe the experimental methods in more detail, with section four presenting our results. The validity and applicability of our work are then discussed.

2. Data

Our subsurface dataset from the Békés Basin consists of Gamma-Ray (GR) logs of both Mezőberény geothermal wells (Figure 2a), along with production logs and permeability from pressure build-up testing in both wells, as well as other data provided by the municipality of Mezőberény. The production logs will be used to estimate the permeability of each sandstone layer in the production intervals of both wells. In addition, three core fragments from the Újfalu Fm. were retrieved from the MOL Plc. core repository: two from the Gyoma-1 well (F3 and F4); one (F1) from the Kond-1 well (Figure 1). GR logs of both these wells are shown in Figure 2b. As is illustrated, both these wells have several cored intervals, typically of ~5 m vertical extent. Core recovery from these wells was very low because of the poor consolidation of the Újfalu Fm, only un-slabbed slices of the cored intervals being available. From the available fragments in the repository, the most sandstone-rich pieces were selected (Figure 3).

During production logging in well VS-1, the total production rate (Q_t) was kept constant at 460 L/min. The log in Figure 2a shows the contributions to the production flow from individual sandstone sections. Twelve production screens have been placed in the production interval of this well, connected by a blind production liner. The log shows that only the ten sandstone layers with lowest GR log readings contribute to the production flow, with individual contribution ranging from 10 to 140 L/min, the combined thickness of all producing layers being 43.3 m. In well T-1 four production screens are installed, of which only three intervals of combined thickness 12.5 m contributed to the production flow. The total production flow was 350 L/min during production logging of this well, the contributions of the three productive layers ranging from 100 to 130 L/min.

Table 1 summarizes these production logging outcomes. The percentage of flow that each layer contributes to the total is presented in the column with header F_i. In addition to production logging, gas-lift production tests were performed with pressure build-up measurements after shut-in periods of ~2 h. This time was too short to determine the skin factor, and only allows rough estimation of the average permeability K_a across the total production interval in each well: 89 mD in well VS-1 and 196 mD in well T-1. Utilizing these values, we estimated the permeability of each production layer assuming a modified version of the standard Theis [30] steady state radial flow equation. Thus:

Q_i = F_i Q_t(1)

where Q_i is the flow rate from layer i, and F_i is the fraction of the total flow rate (Q_t) from that layer, and:

(2)${\mathrm{K}}_{\mathrm{i}}={\mathrm{F}}_{\mathrm{i}}{\mathrm{K}}_{\mathrm{a}}\frac{{\mathrm{h}}_{\mathrm{i}}}{{\mathrm{h}}_{\mathrm{t}}}$

where K_i and h_i are the permeability and thickness of layer i, and h_t is combined thickness of all the producing layers (Table 1).

3. Experimental Rock Characterization Methods

Table 2 summarizes the rock characterization experiments performed on each of the core fragments and the properties thereby determined. The experimental methods are explained in more detail in the following sub-sections.

3.1. X-CT Imaging

X-CT scanning was performed using a Nikon XT H 320/225 system (Nikon Corporation, Tokyo, Japan), equipped with a 225 kV reflection gun and a 2000 × 2000 pixel flat panel photodetector (cell size 0.2 mm × 0.2 mm). The scanning conditions used are summarized in Table 3. Initially, an accelerating voltage of 150 kV was used for each sample with 40 µA current, apart from piece A of sample F4 where a 53 µA current was used. For the first round of scans, the X-ray source to sample distance was set to achieve a minimal voxel size of ~7 µm, as shown in Table 3. A final scan was conducted on a small sub-volume of sample F3 at a higher resolution, after adjusting the accelerating voltage and current (Table 3). No filters were used for these scans. The exposure time for each projection was 1.41 s, repeated for 3141 projections, the total duration for each analysis thus being roughly an hour and a quarter. Figure 4 summarizes the workflow for this imaging and the subsequent image processing.

To ensure accurate quantitative results, high quality images that avoid sources of error, such as artefacts due to the reconstruction process, are necessary. Three-dimensional (3D) volumes were reconstructed from projections using CT Pro 3D software (Nikon Metrology Europe NV, Leuven, Belgium), using an automatic reconstruction tool to find the centre of rotation of each scan and to apply a beam hardening correction [32]. All volumes were reconstructed in 16 bits, giving 65,536 (i.e., 2¹⁶) greyscale values.

3.2. X-CT Image Processing

The 3D image volume was processed using Dragonfly v. 4.1.0.647 software (Object Research Systems Inc., Montreal, Canada) to reconstruct the internal surfaces of the sandstone samples (Figure 4 and Figure 5A). To reduce the effects of ring artefacts and beam hardening at the edges of each image, sub-volumes were created using the “crop” tool (Figure 5B). This mitigated the aforementioned effects and limited the chance of segmenting external air, focusing on the internal pore space. Mineral and pore volumes were separated manually using simple greyscale thresholding, by segmentation of the volumes corresponding to the relative density range of each mineral phase or pore space, a sensitivity analysis being conducted to determine the optimum greyscale threshold for pore space.

3.3. Porosity and Pore Size Distribution with Porosimetry

Small cuboidal pieces, with dimensions of ~3 mm × 3 mm × 3 mm, were cut from core fragments F1, F3 and F4 for porosity measurements. Porosity of these pieces was determined using a Poremaster automated mercury intrusion porosimeter (Quantachrome Corporation, Boynton Beach, FL, USA) [33]. This apparatus progressively applies pressure to force mercury into pore-space in a sample [33,34], using theory for capillary flow into cylindrical pores [35] to determine the cumulative size-distribution of pores in the sample.

3.4. Thin Section Analysis and Scanning Electron Microscopy

Thin sections of each of the samples were analyzed using plane and cross polarized light. These were evaluated based on grain size, sorting, roundness, and mineral composition. Clay minerals were characterized in more detail by Scanning Electron Microscopy (SEM), using a Hitachi SU-6600 field emission scanning electron microscope (Hitachi Ltd., Tokyo, Japan).

3.5. Numerical Flow Experiments

Like other recent studies of sandstones and other lithologies (e.g., [36,37,38]), we used a voxel-based Digital Rock Physics approach for calculation of permeability of sample F3, applying this technique to sub-volumes of the segmented X-CT scan (Figure 4). The X-CT image of this sample was clipped into sub-volume cubes with 1.43 mm (200 voxels) to a side (Figure 5B), 84 sub-volumes in total being modelled. Permeability was calculated independently for flow parallel to each of the x, y, and z axes. This calculation was repeated for three different segmentation thresholds applied to the greyscale X-CT data, resulting in 756 simulations in total. The modelling software used was OpenFOAM v4.1 (The OpenFOAM Foundation Ltd., London, UK). The simpleFoam solver was used; this solves the Navier-Stokes Equations for single-phase incompressible flow. Pore-scale geometry was created using the following workflow:

• Creation of hexahedral meshes for each sub-volume, with every mesh cell directly conformed to a voxel in the X-CT scan. Scan resolution was 7.147 μm, each sub-volume containing 8 million voxels, hence each hexahedral mesh also contained 8 million voxels. All six sides of the bounding box were set as wall-type boundary patches named “Side”.

• Cells that contained a solid voxel were removed from the mesh, leaving only pore space cells. Cell faces exposed by this operation were designated as wall-type boundary patches named “Walls”; these represent the internal boundary between pore space and the rock surface. Total porosity was calculated at this step.

• Inlet and outlet boundary “plates” were added for flow distribution. These plates are analogous to the distribution plates used in rock core holders for experimental measurement of permeability. The boundary patches on these sides were modified to “Inlet” and “Outlet” so that appropriate flow boundaries could be applied. The remaining four external sides of the bounding box stay as no-flow patches named “Side”.

• Regions of pore space not connected to the inlet and outlet plates (i.e., disconnected pore bodies) could not contribute to flow through the sub-volume and were removed. Connected porosity was calculated at this step.

• Steps 3 and 4 were repeated for flow along each axis.

A constant flow velocity of 10⁻¹² m/s into the model domain was specified on the inlet patch, with a fixed pressure of 0 kPa specified at the outlet patch. Walls created when removing solid cells (step 2) were given a no-slip velocity boundary condition as they represent a solid surface at which velocity must be zero. Side walls were given a slip velocity boundary as they represent pore bodies which intersect the sub-volume bounding box. Boundary conditions for each field are listed in Table 4. No turbulence model was included, the inlet flow velocity being set to the specified very low value to ensure flow remained within the laminar regime. The simpleFoam solver was run for each sub-volume, flow direction, and segmentation threshold value until convergence criteria were satisfied. At this point, the average pressure at the inlet was measured and used in Darcy’s equation, along with the specified inlet flow velocity, outlet pressure, fluid viscosity, and dimensions of the sub-volume, to calculate the permeability of that sub-volume.

3.6. X-ray Diffraction

Our XRD analysis aims to characterize sample mineralogy and to identify the presence of authigenic carbonate cement or clay, which may inhibit porosity. Preparation involved gently crushing the samples to powder using a pestle and mortar. XRD patterns were collected using a Rigaku MiniFlex 6G X-ray diffraction apparatus (Rigaku Corporation, Tokyo, Japan) equipped with a D/teX Ultra detector, a 6-position (ASC-6) sample changer and a sealed copper tube (Cu K-α₁ and K-α₂ wavelengths 1.5406 and 1.5444 Å). Diffraction patterns were measured as θ/2θ scans typically over a range of 3° > 2θ > 80°. Data collection and analysis utilized Rigaku SmartLab Studio II software with the Crystallography Open Database [39]. The Reference Intensity Ratio method was used for quantitative analysis of mineralogical composition.

4. Results

4.1. Thin Section and SEM Analyses

Core fragments F1 and F4 can both be classified as siltstones to fine-grained sandstones, and are heavily cemented. Hardly any porosity can be recognized on photomicrographs with plane polarized light (Figure 6A,C). Angular to sub-rounded grains of ~50 µm diameter can be seen; from visual inspection, these consist of ~60% quartz, ~10% potassium feldspar, ~10% micas, and ~20% clay minerals. Linear contact points between the quartz grains and bending within the micas suggest that this sediment has been moderately compacted. The abundant calcite cement, visible with cross polarized light (Figure 6B,D), suggests deposition as a calcite-rich mud, the primary calcite having re-crystallized, following burial, into the calcite cement. SEM analysis indicates the presence of kaolinite, muscovite, and illite clays (Figure 7).

Thin sections indicate that fragment F3 is a medium grained feldspathic quartz arenite sandstone with fine (~125 µm), angular to sub-rounded, well sorted grains (Figure 8A,B). Little to no cementing or clay mineralization is recognized, which explains the unconsolidated nature of this sample. Contacts between grains are dominantly linear and concave, suggesting moderate compaction. SEM images show that the few clay grains present are predominantly kaolinite, small amounts of chlorite, illite, and montmorillonite being also present (Figure 9A–D). A small amount of authigenic quartz was also observed (Figure 9B), which may be a side-product of the transformation of montmorillonite to smectite. This transformation produces silica in solution (e.g., [40]), which could be deposited as authigenic quartz and reduce porosity and permeability in the aquifer. The unconsolidated grains in this sample are most likely easily mobilized and could therefore clog pore space and reduce permeability.

4.2. Mercury Intrusion Porosimetry Measurements

Mercury Intrusion Porosimetry (MIP) measurements were made on one piece of sample F1, three pieces of sample F3, and two pieces of sample F4. As well as porosity, ϕ, density (bulk dry density, ρ_b, including pore space) was also measured as part of this analysis. Grain density was then calculated using the standard formula

(3)${\mathsf{\rho }}_{\mathrm{g}}=\frac{{\mathsf{\rho }}_{\mathrm{b}}}{1-\varphi }$

values thus calculated being included in Table 5. For the sandstone from core fragment F3, the measured porosity varies from 0.296 to 0.355, with a mean value of 0.320 and a standard deviation of 0.025. These values are rather higher than the value of ~0.24 expected, after Almási [23], for sandstone from a depth of ~1900 m. Measured bulk density ranges from 2091 to 2396 kg m⁻³, with a mean value of 2199 kg m⁻³ and standard deviation of 139 kg m⁻³. The measurements on pieces of samples F1 and F4 were combined, as each were siltstones. Porosity varies from 0.063 to 0.135 with mean value of 0.088 and a standard deviation of 0.034. Bulk density varies from 2699 to 3081 kg m⁻³, with a mean value of 2956 kg m⁻³ and a standard deviation of 182 kg m⁻³.

4.3. X-ray Diffraction

The XRD analysis of samples F1 and F3 identified the main mineral phases present: detrital quartz, calcite, and dolomite (Table 6). Groups of minerals such as micas, 1:1 clays, and 2:1 clays were also identified; however, the individual minerals in these groups could not be distinguished by this analysis. Possible candidate minerals are biotite, muscovite or phlogopite for the mica group, kaolinite for the 1:1 clay group, and chlorite, illite, or smectite for the 2:1 clay group. These results indicate that detrital quartz was the dominant mineral type in both samples, accounting for 46.7% of sample F1 and 43.4% of sample F3. Similar proportions of feldspar and micas were identified in each sample, along with pyroxenes and rutile. The dolomite, calcite, and 1:1 and 2:1 clays, identified by the XRD analysis, are diagenetic minerals. These results illustrate that these minerals were present in similar quantities in each sample, the predominant diagenetic mineral being dolomite cement.

4.4. X-CT Results

X-CT scans were run on all three samples, as detailed in Table 3. Except for sample F3 (scan H3), the porosity turned out to be very low, with pore size below the scan resolution. Furthermore, the ‘high resolution’ scan of sample F3 (scan H3-HR) resulted in an indistinguishable porosity distribution compared with the preceding scan H3 of the same sample, indicating that the latter had accurately captured the porosity structure of this sample. X-CT results are therefore only presented here for sample F3 and for scan H3. Figure 10A,B illustrates the pore space segmentation of X-CT scan H3 of sample F3, based on initial visual interpretation. Porosity was calculated for slices of this sample in the x-y plane (i.e., in the horizontal plane, parallel to the bedding) and varied from 0.17 to 0.20 (Figure 10C). The mean porosity for our greyscale cut-off choice is 0.18, its standard deviation being 0.006. Layers of thickness 2–3 mm with different porosity, presumably due to variations in grainsize sorting, can be recognized, porosity being related to grain-size sorting as highlighted by circles 1 and 2 in Figure 10B. In regions with lower porosity (such as circle 1), a larger range in grain size is recognized, whereby the smaller grains reduce the pore space. This ensemble of porosity values is rather less than the ~0.24 expected, after Almási [23], for sandstone from a depth of ~1900 m.

4.5. Numerical Flow Experiments

Numerical flow experiments were carried out on X-CT image H3, derived from sample F3, using the workflow detailed in Section 3.5 (Figure 4). Figure 11 illustrates a representative set of results for the synthetic pressure variations for flow through sub-volumes of this image, in the x-, y-, and z- directions. Three values of greyscale cutoff were adopted, low, mid-range, and high values, which were 7400, 7700, and 8000 for the image as depicted in Figure 5 (and as archived online), representing the range of uncertainty in the identification of the distinction between pore space and mineral grains. Figure 12 illustrates the effect of changing the greyscale threshold, to delineate pore-space, on flow prediction. Many such sub-volumes were analyzed, resulting in a suite of predictions of porosity and permeability, illustrated in Figure 13.

Plotting total porosity against permeability (Figure 13) indicates that these parameters are correlated. It is also apparent that permeability is uniformly highest for flow parallel to the z-axis, and generally slightly higher parallel to the y-axis than parallel to the x-axis. Such variations are evident from the predicted pressure variations being highest in the x-direction and lowest in the z-direction in the images in Figure 11. Figure 12 also indicates that as the greyscale threshold value for segmentation of pore-space increases, the predicted permeability also increases. This is because porosity and pore connectivity increase, creating a greater capacity for flow through the pore network. Figure 14 shows how the specified small increase in the greyscale threshold value increases the connected porosity, resulting in lower flow velocities and fewer bottlenecks, whereas Figure 12 shows how this increase translates into a lower pressure drop across the sub-volume and, hence, higher permeability.

The interdependence between porosity ϕ_i and permeability K_i in a given sub-volume of a granular medium can be modelled [41,42] as a Kozeny-Carman style relationship with percolation threshold ϕ_c (modified from Mavko and Nur [43])

(4)${\mathrm{K}}_{\mathrm{i}}={\mathrm{K}}_{\mathrm{O}}\frac{{({\varphi }_{\mathrm{i}}-{\varphi }_c)}^3}{{(1+{\varphi }_c-{\varphi }_i)}^2}$

where K_O is a reference permeability, used as a fitting parameter. This allows representative values of permeability to be estimated as a function of porosity. The individual estimates of ϕ_i and K_i for different combinations of greyscale cutoff thresholds and flow directions are plotted in Figure 13d for a logarithmic permeability scale, and in Figure 15 for a linear scale, the parameters for fitting the data for each of the combinations of flow directions being listed in Table 7.

To validate this approach, we carried out an equivalent set of numerical flow experiments on a sample of Early Cretaceous Bentheim Sandstone [44] from NW Germany, following the same workflow as in Figure 4 except based on an X-CT image created by Andrew et al. [45]. This approach resulted in a permeability almost identical to their 1900 mD value, experimentally determined in the laboratory at Imperial College, London.

5. Discussion

5.1. Comparison of Results

The permeability values that we have derived from the production tests of wells T-1 and VS-1 at Mezőberény, 18 to 278 mD (Table 1), are low compared to regional averages published in 1994 by Spencer et al. [11], ~70 to ~400 mD (Figure 16). Nonetheless, because Spencer et al. [11] only published average permeability values, it is unclear if our dataset is statistically different from theirs. We have no access to the raw data underlying their study but suggest possible reasons for this difference in results. One possible reason is because of skin formation. However, as already noted, the shut-in periods in between production tests were only two hours long and therefore unsuitable to identify whether skin formation affected the pressure build-up during shut-in. The same short shut-in times also create uncertainty in the determination of the pressure derivative from the pressure build-up curve, which has been used for calculating the average permeability of the production interval. A third possible reason is uncertainty in well completion: if the screened sections of the Mezőberény wells are not accurately aligned with the sandstone intervals then the calculations in Table 1 will be inaccurate. Given these issues, the permeability values calculated in Table 1 must be regarded as only rough estimates.

Our numerical flow simulations have indicated a range of permeabilities of ~150 to ~800 mD in the horizontal direction and ~250 to ~1100 mD in the vertical direction (Figure 13). For comparison, Spencer et al. [11] reported ~100 to ~400 mD in the horizontal direction and ~70 to ~300 mD in the vertical direction (Figure 16). It is thus evident that Spencer et al. [11] regarded this formation as more permeable in the horizontal direction whereas our results indicate greater vertical permeability. However, for both sets of measurements, the difference between horizontal and vertical permeability is small compared with the range of values measured in either direction. Likewise, as already noted, the Almási [23] dataset indicates that the horizontal permeability typically exceeds the vertical permeability by ~50%, both values being typically in the range ~100–200 mD. Furthermore, our estimates of porosity for the high-value greyscale cutoff (Figure 13c,d and Figure 15) are in good agreement with the ~0.24 expected value for a sample from ~1900 m depth, after Almási [23]. Our numerical simulations are based on a single core fragment, F3, from well Gyoma-1 (Figure 1), and cannot necessarily be considered representative of the Újfalu Fm. either at Mezőberény or for the Békés Basin in general. The apparent higher permeability perpendicular to the bedding might be a real feature of the core from the Gyoma-1 well, conceivably through diagenesis, but might potentially instead have arisen from the image processing workflow in which the initial suite of two-dimensional X-ray projections has been reconstructed into a three-dimensional volume, then transformed into a series of slices in planes perpendicular to the z-axis (Figure 4); the Spencer et al. [11] and Almási [23] datasets favour the latter interpretation. Furthermore, the wider range in permeabilities from the flow modelling could be because of the smaller volumes analyzed compared to Spencer et al. [11], which makes it possible to capture high porosity and low porosity sub-volumes separately from each other, whereas in a core-scale experiment these would average out. Indeed, if we were to restrict the permeability analysis to the Kozeny-Carmen fits through the core-averaged porosity of 0.18 (Figure 15), the resulting permeabilities are much more similar to those determined by Spencer et al. [11].

The differences between the permeabilities from the Mezőberény well tests, our numerical flow simulations, and the results published by Spencer et al. [11] might also have a sedimentological origin. Variations in permeability might indeed be anticipated on a variety of scales. First, permeability within fluvial sand bodies has been shown by many workers to vary on a very local scale in different fluvial sedimentary structures (e.g., [46,47,48]). Porosity and permeability can also be expected to differ between turbiditic, deltaic, and fluvial sandstones, so the overall average values presented by Almási [23] do not necessarily represent the properties of any of these sedimentary facies. Second, hydrocarbon fields in the Pannonian Basin are typically found in anticlinal traps or stratigraphic traps. The anticlines are often palaeo-highs, the structural traps being mainly found towards the margins of the Békés Basin [10,11]. In contrast, Mezőberény is located near the centre of this basin (Figure 1). Palaeo-topography may well have impacted sediment supply and the development of accommodation space, which could have affected sedimentation and hence reservoir properties [13]. Suitably-designed production tests would have to be conducted to determine any variation of reservoir properties across this basin. Our results nevertheless raise the possibility that permeability may well have significant lateral variations across the basin.

Comparison is also possible for sandstone sample F3 between the porosities determined from MIP measurements and those deduced from X-CT analysis. The former measurements yielded ~0.32 ± 0.03 (Section 4.2), the latter ~0.180 ± 0.006 (Section 4.4). However, as is illustrated in Table 5, the high porosity from the MIP analysis yields values for the grain density that greatly exceed the ~2650 kg m⁻³ value expected for silica. On the other hand, the mean bulk density of pieces of sample F3, 2199 ± 139 kg m⁻³, adjusts, using our porosity from X-CT analysis and Equation (3), to a grain density of 2680 ± 190 kg m⁻³, a plausible value for silica. We thus conclude that the porosities determined by MIP measurement are unreliable and those determined from X-CT analysis are reliable, justifying our use of these in our investigation of permeability. In this instance, the MIP analysis has not produced reliable results, possibly because the apparatus operated on the basis of theory for cylindrical pores [35], the actual pore spaces being not cylindrical (Figure 10 and Figure 14).

5.2. Injectivity and Productivity Decline

A few months after the Mezőberény wells were completed and the production testing was carried out, the productivity of the injection well (well VS-1) had declined from ~400 L/min to ~40 L/min; as already noted, injectivity decline is common in Hungary. Our thin section analyses support previous studies that suggest scaling and erosion as common causes of these injection problems [19]. Our XRD analysis shows that much of the cement in the analyzed samples is calcite, which could easily dissolve if injection water chemistry were not carefully managed, potentially leading to reservoir erosion and fines migration. In addition, the dissolved calcite could be re-deposited elsewhere, clogging pore space. Another possible cause of the productivity decline is that hydraulic connectivity between both wells is poor [18]. If the sandstone bodies intersected by the wells do not form continuous flow paths between the wells, reinjected water cannot easily flow into and through the reservoir, which increases the pressure required to sustain flow rates [49]. Connectivity issues are common in sedimentary reservoirs, especially if the sandstone content is low as the GR logs of well VS-1 and T-1 indicate (Figure 2) [50,51]. Ainsworth et al. [52] reported an example of the impact of poor reservoir connectivity on recovery efficiency in the Sirikit hydrocarbon field in Thailand, which may be an analogue of the Újfalu Fm.

Another injectivity-reducing factor was recently recognized in the Mezőberény wells; in 2017, biofilm was found in the reinjection well VS-1 following maintenance work. We consider it unlikely that this biofilm caused the early injectivity decline in this well, but we mention it here to give a full overview of injection issues in the Újfalu Fm. Microbiology might both reduce injectivity as well as enhance scaling and reservoir erosion because it could influence pH and water chemistry. The microbial activity at Mezőberény might have been introduced by circulating insufficiently cleaned surface water during maintenance; alternatively, autochthonous microbes might have been activated by injecting water that was chemically different from the original formation water [53]. Indeed, Osvald et al. [54] have presented a detailed analysis of the biofilm-related issues that have affected the geothermal energy project at Hódmezövásárhely in SE Hungary (~25 km NE of Szeged; Figure 1).

5.3. Application of Our Results

Our experiments provide detailed information on the mineralogy, porosity, and permeability of the analyzed sandstone core fragments. However, upscaling these results for reservoir modelling purposes remains a challenge because only three core fragments have been analyzed. To generate actual distributions of hydraulic properties of Újfalu Fm. sandstones, more experiments like those described will be required on core fragments from other sandstone intervals and also from different wells. Our experiments nonetheless indicate how the understanding of these properties could be improved by future work.

Another difficulty limiting the application of our work in future reservoir modelling or feasibility studies is that the properties that we have measured cannot be simply related to gamma ray log signals, such as those depicted in Figure 2. This is because the core fragments from each of the cored intervals were stored in bags, without keeping track of the original stratigraphic order. Furthermore, for core interval 1 of well Kond-1, no GR log is available at all (Figure 2). An analysis, using the methods in our present study, should be carried out on a suite of core samples that are accurately located within a well for which a GR log is also available. It should ultimately be feasible to convert the general association, evident in Figure 2, between low GR readings and favourable reservoir properties, into a quantitative calibration tool where GR readings can act as a quantitative proxy for reservoir porosity and permeability.

In principle, exploration techniques developed for the petroleum industry are transferable to the geothermal sector. However, the lower value of heat compared with hydrocarbons precludes use of expensive exploration methods (e.g., [55]). X-CT analysis followed by numerical flow experiments is not yet recognized as a standard approach to determination of reservoir hydraulic properties [56,57], being considered an emerging technology (e.g., [37]), facilitated by the growing availability of both X-CT scanners in universities and other laboratories [58] and high memory computers, capable of the necessary processing. Although this approach has already been applied to the analysis of sandstone samples (e.g., [36,37]), we are not aware of its previous use in the geothermal sector. Many geothermal projects (e.g., [59]) indeed depend on existing reservoir permeability data from studies of core on behalf of the petroleum sector; where new data are needed, the options have hitherto been limited. In one recent case study [60], where (as at Mezőberény) existing core was scarce, only a few laboratory measurements of permeability were possible; much more information could have been obtained using X-CT. In another [61], a new borehole was cored to provide core for petrophysical analysis; with use of X-CT, maybe on nearby outcrop samples, the additional cost of drilling with coring, over drilling without coring, could have been avoided. X-CT analysis followed by numerical flow experiments indeed warrants recognition as a potential contributor to reducing the cost of geothermal exploration, enhancing the economic viability of projects, and thus facilitating decarbonization of heat supply.

6. Conclusions

Well tests and production logging in geothermal wells at Mezőberény in SE Hungary indicate that the initial permeability of sandstone layers in the Újfalu Fm. ranged from several tens of mD to ~250 mD, below previously published regional averages. In contrast, our numerical flow simulations have indicated a range of permeabilities of ~150 to ~1100 mD, encompassing and exceeding the previously published values. The production logs show that some sandstone layers, in which the production well was completed for production, proved to have zero productivity, highlighting the value of production logging for detailed reservoir performance analysis, and identifying net reservoir thickness in the Újfalu Fm. X-CT scans and thin section analyses show that Újfalu Fm. sandstone has porosity ~0.180 ± 0.006. This lithology can be poorly consolidated, making it prone to reservoir erosion and creating challenges for well cleaning treatments, as dissolution of clogged material in pore space could easily cause further reservoir erosion. Finer-grained sandstone intervals can be completely cemented with calcite, only intervals with the lowest gamma ray readings proving productive by production logging. By eliminating the need for expensive bespoke coring to obtain reservoir properties, this workflow, including X-CT scanning followed by numerical flow simulation, will contribute to reducing the capital cost of developing geothermal energy projects, thus facilitating decarbonization of global energy supply.

Author Contributions

Conceptualization, H.M. and R.W.; methodology, J.M. and H.M.; software, J.M.; validation, C.J.L.W. and R.W.; formal analysis and investigation, C.C., S.M.W., A.W. and D.W.; resources, J.M.; data curation, C.J.L.W., C.C. and J.M.; writing—original draft preparation, C.J.L.W., C.C. and J.M.; writing—review and editing, R.W.; visualization, C.J.L.W., J.M. and R.W.; supervision, H.M., N.M.B. and R.W.; project administration, N.M.B. and R.W.; funding acquisition, R.W. All authors have read and agreed to the published version of the manuscript.

Funding

This research and the APC were funded by European Commission Horizon 2020 research and innovation programme under grant agreement No. 691728 (DESTRESS).

Data Availability Statement

Existing data used in this study are available from the sources stated in the text. Our X-CT images and porosity and permeability estimates from X-CT analysis followed by numerical flow experiments are available from the University of Glasgow online repository at doi:10.5525/gla.researchdata.1140. The remaining data generated in this study are reported in the body of this article.

Acknowledgments

We are grateful to the municipality of Mezőberény for providing data on their geothermal wells, as part of the DESTRESS Horizon 2020 project, to Kirsten Nowak for collecting these data, and to Abel Marko and Máté Osvald for translating them from Hungarian. We are also grateful to MOL Plc. (Magyar OLaj- és Gázipari Részvénytársaság, Hungarian Oil and Gas Public Limited Company) for access to their core repository. We thank the staff of the Advanced Materials Research Laboratory of the University of Strathclyde for their assistance with the experimental work: Tiziana Marocco and Fiona Sillars for the hydraulic property analysis; and Alice Macente for the X-CT data acquisition and reconstruction. The Scottish Oil and Gas Technology Centre (formerly the Oil and Gas Innovation Centre) is also acknowledged for funding the capital cost of this X-CT equipment. Matthew Andrew kindly provided X-CT imagery of the Bentheim Sandstone, which was used to validate our numerical flow experiments. We also thank Object Research Systems for granting an academic license for their Dragonfly image analysis software and Panterra Geoconsultants for granting an academic license for their Cyclolog ENRES INPEFA software. We are also grateful to the three anonymous reviewers for their thoughtful and constructive comments.

Conflicts of Interest

The authors declare that no conflict of interest exists. The funder had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Location of the study area. (A) Stratigraphic cross-section across the Békés Basin and neighbouring Mako Trough, sub-basins of the wider Pannonian Basin. (B) Map of Békés province in SE Hungary, showing the study area, with hydrocarbon and geothermal wells and seismic lines. (C) More detailed map of the study area, showing the Mezőberény well doublet, hydrocarbon wells that have yielded core fragments, and seismic lines. Inset shows locations of parts (A–C) in SE Hungary. Modified after Figure 1 of Brehme et al. [15].

Figure 2. Gamma Ray logs of Békés Basin wells. (a) Mezőberény geothermal production well T-1 (well B-115) and re-injection well VS-1 (well K-116), with production data also shown. (b) Hydrocarbon wells Kond-1 and Gyoma-1.

Figure 3. Descriptions of recovered core fragments.

Figure 4. Summary of the workflow for X-CT imaging, followed by image processing and numerical flow simulation. Modified after Figure 1 of [31].

Figure 5. X-CT imagery. (A) Section of X-CT scan segmented to show solids as grey and pore space as black. (B) Sub-volume cube.

Figure 6. Photomicrographs of samples F1 and F4. (A) Sample F1 under plane polarized light. (B) Sample F1 under cross-polarized light. (C) Sample F4 under plane polarized light. (D) Sample F4 under cross-polarized light.

Figure 7. SEM images of sample F4. (A) Wide, low-resolution view of part of this sample, showing diverse secondary mineralization. (B) Higher resolution image of a smaller part of the sample, showing secondary kaolinite.

Figure 8. Photomicrographs of sample F3 under (A) plane polarized and (B) cross-polarized light.

Figure 9. SEM images of sample F3. (A) Showing diverse secondary minerals. (B) Showing authigenic quartz. (C) Showing quartz and kaolinite. (D) Showing kaolinite.

Figure 10. Analysis of X-CT results for sample F3 based on analysis H3. (A) X-y plane of X-CT scan. (B) The image in (A), labelled to show segmented pore space, also showing variations in grain size sorting, indicating the range of porosity values and their vertical trend. (C) Calculated porosity for multiple slices through this sample parallel to the x-y plane, showing the variability in porosity that contributes to the determination of the mean value.

Figure 11. Simulated pressure variations across a sub-volume of sample F3, calculated for flow parallel to the x, y, and z axes.

Figure 12. Simulated pressure variations across a sub-volume of sample F3, calculated for flow parallel to the z axis, taking different grey-scale thresholds as delineating pore-space from filled space. (A) Using the mid-range grey-scale threshold. (B) Using the high grey-scale threshold, which predicts more pore space enabling flow to occur with much smaller pressure variations, giving a higher permeability.

Figure 13. Total porosity and permeability for each analyzed sub-volume, for simulated flow parallel to the x, y and z axes, for the three greyscale cutoffs for pore space. Thus, for example, ‘mid-Z’ means mid range cutoff and flow in the z direction. (a) Low cutoff, 7400. (b) Mid range cutoff, 7700. (c) High cutoff, 8000. (d) The combined datasets, along with Kozeny-Carman type curve fits, after Equation (4), using parameters in Table 7, for flow in the z-direction (i.e., vertical), the x and y directions, and in all three directions. The data presented in this Figure are available at doi:10.5525/gla.researchdata.1140.

Figure 14. Simulated velocity distributions within a sub-volume of sample F3, again calculated for flow parallel to the z-axis, taking different greyscale thresholds as delineating pore-space from filled space. Colourless areas denote volume interpreted as rock; coloured areas denote flow velocities in volume interpreted as pore space. (A) Using the mid-range greyscale threshold. (B) Using the high greyscale threshold, which predicts opening of more pore space enabling flow to occur much more freely through the model sample.

Figure 15. Enlargement of Figure 13d with linear scale, showing estimates of permeability and porosity for different sub-volumes of sample F3, greyscale cutoffs for pore space, and modelled flow directions, fitted using modified Kozeny-Carman curves of the form of Equation (4).

Figure 16. Comparison of results of this study with those of Spencer et al. [11].

Table 1

Production logging of wells VS-1 and T-1.

Production	Depth	Depth	Thickness	Qi	Fi	Ki
Interval, i	Top [m]	Base [m]	Hi [m]	[L/min]	[%]	[mD]
Well VS-1
1	1643.9	1646.5	2.6	20	4.4	65
3	1665.8	1669.7	3.9	20	4.4	44
4	1745.0	1755.0	10.0	110	23.9	92
5	1793.0	1794.7	1.7	30	6.5	147
6	1810.7	1817.9	7.2	15	3.3	18
7	1856.0	1857.3	1.3	20	4.3	128
9	1882.0	1886.5	4.5	25	5.4	46
10	1891.0	1893.7	2.7	70	15.2	217
11	1898.0	1899.8	1.8	10	2.2	47
12	1920.0	1927.6	7.6	140	30.4	154
Total			43.3	460
Well T-1
1	1828.5	1833	4.5	130	37	201
2	1907.5	1910.5	3.0	120	34	278
3	1942.5	1947.5	5.0	100	29	142
Total:			12.5	350

This analysis is based on Equations (1) and (2). K_i is the estimated permeability of the ith production interval.

Table 2

Summary of experimental analyses.

Sample	F1	F3	F4
Facies	Siltstone	Med. Grained Sst	Fine Grained Sst.
Analysis				Property
X-CT imaging	Y	Y	Y	Porosity
X-ray Diffraction (XRD)	Y	Y		Minerology
Thin section	Y	Y	Y	Minerology
Scanning Electron Microscope (SEM)	Y	Y	Y	Minerology
Mercury Intrusion Porosimetry (MIP)	Y	Y	Y	Porosity, density
Numerical flow simulation		Y		Permeability

Table 3

X-CT scan conditions.

Sample	Notes	Scan	Voxel Size (µm)	Scan Size (Voxels)	Voltage (kV)	Current (µA)
F1		H1	7.47888	1684 × 1436 × 1592	150	40
F4	(1)	H2	8.78948	1494 × 1556 × 1779	150	53
F3		H3	7.14653	1590 × 1751 × 1866	150	40
F4	(2)	H4	6.52282	1609 × 1370 × 1936	150	40
F3	(3)	H3-HR	2.78965	2000 × 2000 × 1399	140	19

Notes: (1) piece ‘A’ of sample F4; (2) piece ‘B’ of sample F4; (3) ‘High resolution’ (HR) analysis of a sub-volume of sample F3.

Table 4

Boundary conditions in numerical flow experiments.

Boundary Name	Type	Velocity Field	Pressure Field
Inlet	patch	fixed velocity	zero gradient
Outlet	patch	zero gradient	fixed pressure
Side	wall	slip	zero gradient
Walls	wall	no-slip	zero gradient

Table 5

Mercury intrusion porosimetry measurements.

Sample	Piece	Porosity	Bulk Density	Grain Density
			ρb (kg m−3)	ρg (kg m−3)
Sandstone
F3	1	0.296	2396	3403
F3	2	0.355	2091	3242
F3	3	0.310	2111	3059
Siltstone
F1	1	0.135	2699	3120
F4	1	0.066	3088	3306
F4	2	0.063	3081	3288

Table 6

Mineralogical compositions from XRD.

Sample	Mineralogical Composition (%)
	Quartz	Calcite	Dolomite	Mica	Feldspar	1:1 Clay	2:1 Clay	Rutile	Pyroxene
F1	46.7	4.7	29.5	2.4	10.0	1.0	3.5	1.6	0.3
F3	43.4	5.4	25.0	2.8	10.4	1.0	3.3	0.7	8.0

Table 7

Kozeny-Carman parameters for sample F3.

Flow Direction	KO (10−12 m2)	KO (mD)	ϕ c
z	104.36	105,743	0.0287
x, y and z	71.309	72,254	0.0208
x and y	77.006	78,026	0.0329

1 mD ≡ 9.869233 × 10⁻¹⁶ m².