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

: 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. How-ever, 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.


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 −2 ) and geothermal gradient (~45 °C km −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.

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 (Qt) 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 Fi. 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 Ka 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: where Qi is the flow rate from layer i, and Fi is the fraction of the total flow rate (Qt) from that layer, and: where Ki and hi are the permeability and thickness of layer i, and ht is combined thickness of all the producing layers (Table 1). This analysis is based on Equations (1) and (2). Ki is the estimated permeability of the i'th production interval. 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.

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 subvolume 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 16 ) greyscale values.

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 (Figures 4 and 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.

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, Florida) [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.

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).

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 Open-FOAM 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: 1. 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". 2. 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. 3. 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". 4. 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. 5. Steps 3 and 4 were repeated for flow along each axis.
A constant flow velocity of 10 −12 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.

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-α1 and K-α2 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.

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.

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 b g ρ ρ = 1-φ (3) 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 −3 , with a mean value of 2199 kg m −3 and standard deviation of 139 kg m −3 . 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 −3 , with a mean value of 2956 kg m −3 and a standard deviation of 182 kg m -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.

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.

Numerical Flow Experiments
Numerical flow experiments were carried out on X-CT image H3, derived from sample F3, using the workflow detailed in subsection 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 xaxis. 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 Ki 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]) where KO 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 Ki 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.

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 (Figures 13c,d and 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 −3 value expected for silica. On the other hand, the mean bulk density of pieces of sample F3, 2199 ± 139 kg m −3 , adjusts, using our porosity from X-CT analysis and Equation (3), to a grain density of 2680 ± 190 kg m −3 , 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 (Figures 10 and  14).

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).

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.

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.