Experimental Rock Characterisation of Upper Pannonian Sandstones from Szentes Geothermal Field, Hungary

Abstract

The Upper Pannonian (UP) sandstone formation has been utilised for thermal water production in Hungary for several decades. Although sustainable utilisation requires the reinjection of cooled geothermal brine into the host rock, only a fraction of the used water is reinjected in the country. UP sandstone formation is reported to exhibit low injectivity, making reinjection challenging, and its petrophysical properties are poorly known, which increases uncertainty in designing operational parameters. The goal of the study is to provide experimental data and to gain a better understanding of formation characteristics that control injectivity and productivity issues in Upper Pannonian sandstone layers. Petrographical characterisation and petrophysical laboratory experiments are conducted on cores retrieved from two wells drilled in the framework of an R&D project at the depth of between 1750 m and 2000 m. The experiments, such as grain density, porosity, permeability, and ultrasonic velocity, as well as thin section, grain size distribution, XRD, and SEM analyses, are used to determine Petrophysical Rock Types (PRT) that share distinct hydraulic (flow zone indicator, FZI) and petrophysical characteristics. These are used to identify well intervals with lower potential for injectivity issues. The results imply that fines migration due to formation erosion is one of the key processes that must be better understood and controlled in order to mitigate injectivity issues at the study area. Future investigation should include numerical and experimental characterisation of formation damage, including water–rock interaction tests, critical flow velocity measurements, and fines migration analysis under reservoir conditions.

1. Introduction

The Upper Pannonian (UP) sandstone formation in Hungary has been utilised for thermal water production, especially in the Szentes Geothermal Field, for over 60 years [1]. Although sustainable utilisation requires the reinjection of cooled geothermal brine into the host rock, less than 10% of all geothermal wells have been used as reinjection wells in the country [1]. This is linked to economic constraints posed by high-pressure injection technology employed by the oil and gas industry, as well as several unsuccessful reinjection operations in Pannonian sandstone formations [2,3,4].

Several underlying mechanisms were proposed that might cause reinjection issues, but no effective mitigation strategy has been developed yet [5]. Markó et al. [6] have recently proposed a methodology for the systematic identification of potential reasons for the low injectivity of sandstone aquifers with suggested methods for eliminating them. According to this study, possible processes that limit reinjection can occur in the near borehole area, or be related to reservoir or regional scale hydraulics. At borehole scale, local clogging processes such as particle migration, e.g., sand production due to formation damage [7,8,9], and mineral precipitation, e.g., scaling [10,11], as well as microbial activity, e.g., biofilm formation [12,13], are possible. At reservoir scale, inadequate performance may be associated with limited reservoir extent [14] and low permeability as well as reservoir performance [15]. At regional scale, the potential presence of an overpressured zone may be considered [16,17]. The injection problems may arise from a combination of these processes as well. Using the workflow proposed by Markó et al. [6], Brehme et al. [18] report on the successful implementation of this approach in a geothermal well at a depth of approx. 2000 m near Mezőberény in the Békés Basin, South-East Hungary. They conclude that low injectivity of Újfalu Formation sandstone rock is associated with low reservoir permeability and the precipitation of carbonates, iron, and manganese minerals. Injectivity enhancement is achieved by the combination of tailored hydraulic and chemical stimulation.

Besides investigating possible sources for reduced injectivity, a major challenge in planning field operations is the limited availability of petrophysical properties for UP sandstone formations [18,19,20,21]. Willems et al. [20] have recently proposed a methodology for filling this data gap based on laboratory tests conducted on Újfalu Formation core fragments from two legacy wells in the Békés Basin and subsequent numerical flow simulation. The experiments include petrographic methods, such as thin-section and X-ray Diffraction (XRD) analysis, and petrophysical techniques including helium gas porosimetry and X-ray Computed Tomography (X-CT) imaging, as well as numerical flow simulation for determining permeability indirectly. Although the applicability of upscaling these results for reservoir scale studies is limited by the very small number of core fragments analysed, i.e., solely three samples, this is enhanced by numerical flow simulation for determining permeability indirectly.

Despite the previous efforts and proposed workflows discussed above, geothermal doublets drilled in UP sandstone formations still experience injection problems. Therefore, there is a need for better understanding of both formation characteristics and the underlying mechanisms resulting in low injectivity in order to develop effective mitigation strategies. The research and development (R&D) project “Development of a well completion technology for sustainable and cost-effective reinjection of thermal water” aims at the multiscale characterisation of UP sandstone formations and development of a targeted methodology for maintaining or enhancing the injectivity of geothermal wells with a focus on the Szentes Geothermal Field in Hungary [18]. The project includes a series of laboratory experiments, the drilling of new wells targeted at these formations including geophysical logging and hydraulic investigations, and the recompletion of an existing well for reinjection for Frac&Pack stimulation. We note that the construction of the well SZT-1 and the recompletion technology for the existing well SZT-VIII (also referred to as K-666) are presented in Farkas et al. 2022 [15] in detail.

In this paper, we present the results of the laboratory experiments UP sandstone samples in the framework of the R&D project. In this study, we focus on the core scale and microscale features of the investigated samples in order to give an insight into the petrophysical and petrological characteristics of UP sandstone rocks. The laboratory investigation includes petrographic characterisation, such as thin section, grain size distribution, XRD and scanning electron microscope (SEM) analyses, and petrophysical experiments, i.e., grain density, helium gas porosity, and permeability, as well as ultrasonic wave velocity measurements. The possibility of conducting these tests on several cores retrieved from newly drilled exploration wells allows extending the limited public dataset on mineralogical and petrophysical characteristics of UP sandstones. Furthermore, the analyses may contribute to the understanding of characteristics that control decline in injectivity or productivity in UP unconsolidated sandstone reservoir. We apply the Petrophysical Rock Typing (PRT) technique, which allows classifying rocks that share similar hydraulic and petrophysical properties for identifying well intervals with lower potential for injectivity problems.

This paper is structured as follows: in Section 2, the Szentes Geothermal Field is presented. In Section 3, the experimental methods are described. In Section 4, the petrographic and petrophysical results are presented. In the same section, the results are compared; their applicability and suggestions for future work are also discussed. In Section 5, conclusions are drawn.

2. Szentes Geothermal Field

2.1. Field History

The Szentes Geothermal Field is located in South-East Hungary on the left bank of the Tisza river (Figure 1). It is one of the most intensively utilised geothermal areas in Hungary, with 40 active wells producing more than 5.5 million m³ of hot water per year [1]. The produced thermal water is utilised for district heating and balneology, as well as for agricultural purposes. All the produced water is discharged in surface water [3].

The latest study on the production history of the Szentes Geothermal Field is published by Bálint and Szanyi 2015 [1], who provide an in-depth overview of field development and hydraulic characteristics, such as production history of the wells based on previous hydraulic test reports [22,23,24] and the latest hydraulic test campaign, conducted in 20 wells between 2009 and 2010 [2]. They point out that continuous production over decades without reinjection results in a significant drop in production rate, i.e., approx. 7.6 million m³/year in 1971–1972 vs. 5.5 million m³/year in 2009–2010, and a pressure drop of 1.5 to 4 bar with respect to hydrostatic pressure, both factors contributing to the decline in injectivity of the doublet.

2.2. Geological Setting

The area is part of the northern wedging of Makó Trough, where the depth of the Pre-Neogene basement ranges in depth from 3000 to 5000 m [1,25]. The generalised chrono- and lithostratigraphy of the study area is shown in Figure 2 and the depth map of the basement top is illustrated in Figure 3. Basin subsidence began in the Miocene, with the highest rate in the Pannonian (Late Miocene to Pliocene), resulting in sediment deposition with a thickness of more than 4000 m. Lower Pannonian (LP) sandstone formations, i.e., Endrőd Fm., Szolnok Fm., and Algyő Fm., are characterised by clay–marl and very fine-grained powdered quartz layers. The lower part of the Upper Pannonian (UP) sandstone formations is characterised as sandstone with clay–aleurite streaks or laminae as a result of deposition in delta plain, moor, and smaller bay environments. The upper 300–400 m of these sandstone layers are described as loose, poorly consolidated sandstone. The top of the Pannonian sandstone formations in the study area is between 2000 and 2500 m, where lower elevations are located towards Szegvár, south from Szentes (Figure 4). The UP sandstone layers are covered by Pliocene and Pleistocene sediments (Figure 2). We note that, according to the latest stratigraphic nomenclature, the Zagyva and Újfalu Fms. are referred to as Transdanubian Formation Group, and the Algyő, Szolnok, and Endrőd Fms., as well as Tótkomlós marl, are referred to as Alföld Formation Group [26]. For the sake of clarity, we use the term UP sandstone formation as lithofacies associations of Zagyva and Újfalu Fms. in this paper.

In most of the geothermal wells at the Szentes Geothermal Field, the production intervals are perforated in the Újfalu Fm. (Section 2.3); thus, we focus on this rock formation. The formation is composed of fine and medium sandstone intercalated by thin marl and siltstone layers. Sandstone bodies are of estuary bar, delta branch riverbed filling, and crevasse splay origin. The intercalating layers are associated with oxbow lake environment [25].

2.3. Geothermal Reservoir Characterisation

Based on production history and well test analysis of the 40 active wells in the Szentes Geothermal Field, three aquifer layer groups can be defined [1]. Most of the wells have a completion with production intervals in the Újfalu Fm. A stratigraphic cross-section across the study area with wells and their production intervals is shown in Figure 5. The upper aquifer layer group, level A, consists of wells having a completion with production intervals in the Újfalu and Zagyva Fms. between the depth of 1500 and 1800 m with an average permeability of 1500 mD. The middle aquifer layer group, level B, includes wells with production intervals between the depth of 1800 and 2000 m, mainly in Újfalu and partly in Zagyva Fm. Rock, with an average permeability of 500 mD. The lower aquifer layer group, level C, includes wells below the depth of 2000 m entirely in Újfalu Fm. with an average permeability varying between 1000 and 2000 mD. Thermal water production is dominated by wells screened in level B. A summary of 14 wells including well ID, location coordinates in the Hungarian national projection system (EOV), drilling year, depth, screening, and production rate, as well as bottom-hole temperature, is presented in

Table 1. Characteristics of wells close to “SZT-1” and “SZSZT-IX” located in Szentes (based on Bálint and Szanyi [1]).

Local Well ID	National Well ID	EOV Y (m)	EOV X (m)	Drilling Year	Depth (m)	Production Interval (m)	Bottom-Hole Temperature (°C)	Production Rate (m3/year)
SZT-I	K-498	747,458	149,539	1964	1995	1800–1975	85	154,860
SZT-II	K-562	747,489	149,493	1970	1800	1640–1793	82	150,850
SZT-III	K-563	746,484	148,355	1970	1992	1678–1936	78	170,150
SZT-IV	K-586	747,234	149,283	1972	2303	2060–2235	96	178,604
SZT-V/1	K-640	747,855	150,857	1979	2240	2040–2210	94	134,800
SZT-V/2	K-641	747,800	150,800	1979	2000	1785–1993	84	155,100
SZT-VI/1	K-642	749,981	150,148	1978	2398	2046–2255	97	172,500
SZT-VI/2	K-643	749,982	151,063	1978	1998	1694–1989	86	189,000
SZT-VII/1	K-644	747,112	152,442	1979	2257	2053–2205	96	156,900
SZT-VII/2	K-639	747,111	152,480	1979	1806	1534–1754	76	174,300
SZT-VII/3	K-645	747,101	152,539	1980	1998	1800–1998	80	194,100
SZT-VIII	K-666	749,166	151,854	1988	2300	2004–2143	90	128,900
AL/1	K-561	741,347	151,764	1969	2050	1801–2019	85	256,000
AL/2	K-578	742,231	151,198	1971	2401	2135–2401	94	256,000

.

3. Sample Collection and Experimental Methods

3.1. Core Sample Collection and Description

In the framework of this R&D project, two vertical exploration wells, referred to as “SZT-1” and “SZSZT-IX”, were drilled in Szentes in 2020 to retrieve core samples for laboratory experiments and to conduct a long-term reinjection test at well SZT-1. The coring intervals were determined based on the stratigraphy of the offset wells, K-564 and K-515, and the seismic interpretation of the study area reported by Bereczki et al. 2020 [29]. Cores were collected between approximately 1740 m and 1970 m depth in order to penetrate level B Újfalu sandstone layers approximated from the nearest offset wells, K-564 and K-515.

Table 2. Well parameters and coring properties of boreholes SZT-1 and SZSZT-IX.

Local Well ID			SZT-1	SZSZT-IX
National Well ID			K-712	K-707
EOV Y (m)			748,464	747,995
EOV X (m)			149,886	149,416
Drilling year			2020	2020
Total depth—MD (m)			2000	2009.2
Bottom-holetemperature (°C)			92.8	88.0
Screen intervals (m)		Top	1934.0	1932.4
		Bottom	1981.4	1997.3
Coring intervals	1	Top	1740 m	1742.5 m
		Bottom	1749.25 m	1755.38 m
	2	Top	1930 m	1835 m
		Bottom	1972.5 m	1840.33 m
Average core recovery			85.2%	85.1%

summarises the main parameters of the wells. Figure 6 shows that the total length of collected samples is equal to 52.17 m with an average core recovery rate of 85 % from the two boreholes.

After core retrieval and documentation (Figure 6), the logged samples were packed in cling film in order to preserve their moisture content. These were transported in wooden boxes to the laboratory.

Light grey, very fine-fine grained, poorly cemented, micaceous sandstone is the most common rock type in the investigated depth intervals. It is carbonate cemented very well in some places. Very fine-fine grained sandstone usually appears above fine-medium grained and medium-coarse grained sandstone and below very fine grained sandstone, with siltstone forming fining upward sequences. Coarsening upward sandstone sequences also appear in the strata. The sandstones consist of quartz, feldspar, carbonate rock debris, mica, and clay minerals (mainly kaolinite and illite), as well as coalified plant fragments. Siltstone, argillaceous marl, and coaly argillaceous marl appear between the sequences, which is a hint toward a low-energy environment. Based on the depositional environments, several sedimentary structures can be observed in the cores, as shown in Figure 7.

Figure 8 illustrates representative identified depositional facies in well SZT-1 showing typical channel–overbank sequence of a meandering channel based on analysis of sedimentary structures of core samples and the shape of the GR well log. For more details of the stratigraphical model of the Upper Pannonian sandstone sequence in the study area, the reader is referred to [30].

3.2. Sample Preparation

Grain density, porosity, permeability, and ultrasonic wave velocity measurements were carried out on cylindrical rock samples. The plugs were drilled with a diameter of 1.5”, both parallel and perpendicular to the core axis. After the drilling procedure, samples were saw-cut and the end faces of the samples were carefully polished with a grinder machine to reach the desired parallelism in accordance with ASTM and ISRM standards [31,32].

The plugs were dried at a temperature of 60 °C to preserve the chemically bound water in the lattice of clay minerals and then were stored in a desiccator between each measurement. In order to prevent sample contamination, coupling media were not used for ultrasonic velocity measurements.

Grain size distribution measurements, thin section analysis, XRD, and SEM measurements were carried out on the remaining rock slabs of the plug samples. To reduce the charging effect in the case of SEM imaging, the test specimens were coated with gold.

3.3. Laboratory Experiment Methods

First, thin sections of 15 samples from borehole SZSZT-IX were analysed. After that, petrophysical measurements, including grain density, porosity, permeability, and ultrasonic measurements, were performed on the 1.5” plug samples. Grain size distribution was measured on the remaining rock slices after slabbing the plug samples.

Based on the petrophysical data processing results, representative samples of Petrophysical Rock Types (Section “Petrophysical Rock Typing”) were selected according to stratified sampling strategy [33] for further petrographic analysis, including SEM and XRD, to investigate textural features.

The applied methods are presented based on their nature in subsequent Section 3.3.1 and Section 3.3.2.

3.3.1. Petrographical Characterisation

The thin sections of the samples were analysed with a Carl Zeiss polarized light microscope using plane-polarized and cross-polarized light. These were evaluated based on grain size, sorting, roundness, and mineral composition.

Grain size distribution was measured by laser diffraction method using Cilas 1180 device. The filter cake was carefully removed from the core surfaces. The samples were disaggregated using distilled water. Prior to measurement, each sample was ultrasonicated for 180 s under stirring conditions and also during the measurement in order to ensure sample dispersion. This measurement was performed at least 3 times on a sample.

The XRD patterns on the sandstone samples were collected using Cu-Kα radiation (40 kV, 15 mA) with a Rigaku MiniFlex 600 (Rigaku, Tokyo, Japan). Scans were made at room temperature from 5 to 70° 2θ, with a step of 0.02/s. XRD scans were evaluated for quantitative phase composition using a full profile fit procedure. The total amount of identified (crystalline) phases is taken as 100%. Due to the unknown proportion of amorphous components, the phase percentages reflect only relative abundances. The measurement uncertainties are ±1%, due to the precise sample preparation and measurement.

SEM imaging was conducted with a Jeol JSM-IT500HR (Jeol, Tokyo, Japan) instrument. Measurements were performed in a high vacuum chamber with a beam voltage of 5.0 kV.

3.3.2. Petrophysical Experiments

Grain Density

Matrix volume was measured by a Quantachrome Pentapyc 5200e (PPY-30T) instrument. This test follows the principle of the Boyle–Mariotte Law. A known amount of He flows through on a given pressure from the reference cell with VR volume to the sample chamber. The volume of the sample chamber (VC) is determined by calibration of the instrument with stainless steel reference spheres at a given temperature before the measurement. Seven measurements were carried out on each sample but the average grain volume was calculated from the last five values. Measurements were performed in a tempered thermostat at a constant temperature of 25 °C. For grain density calculation, the weight of the sample was measured by an analytical balance with 0.1 mg accuracy. The bulk volume and porosity of the plugs were calculated from geometrical data of 3D scanning.

Porosity and Permeability

He gas porosity and permeability under reservoir pressure conditions were measured by Vinci Technologies COREVAL-700 gas permeameter. The plug samples were measured after He pycnometry. The plugs were held in an isostatic core holder during the tests. The applied confining pressure was 210 bar for each sample at lab temperature. The method used for determining porosity in this case is called “Boyle’s Law Single Cell Method for direct void volume measurement” [34]. He gas permeability was based on “Transient pressure technique for gases: Pressure-Falloff, Axial Gas Flow measurements” [34]. This technique has a useful permeability range of 0.001 to 5000 mD. The measured gas permeability was corrected for the Klinkenberg effect to obtain water permeability.

We note that several historical permeability data for similar formation sandstone rock are available in [35,36]. However, these present only a compilation of datasets instead of original experimental data. For more details on these data, we refer to Willems et al. [20].

Ultrasonic Velocity

The ultrasonic velocity of compressional and shear waves was measured by SRL A1000 instruments using a pulse-transmission technique [31,32]. In this case, two transducers were placed on the end faces of the samples. The frequency of the transducers used for measurements was centered around 1 MHz, both for compressional and shear waves. Travel times for velocity data were determined with „first-break” record using a modified Akaike Information Criterion algorithm [37].

Experimental Setup for Porosity, Permeability, and Ultrasonic Velocity Measurements

Samples were put into a high-pressure isostatic core holder for porosity, permeability, and ultrasonic velocity measurements to mimic in situ reservoir pressure conditions. The applied pressure was calculated using the equation for linear poroelasticity [38]:

$$P_{eff}=S_L-\alpha P_p$$

(1)

where P_eff is the effective pressure, S_L is the uniform lithostatic stress, α is Biot’s coefficent, and P_p is pore pressure. Since α is not known for UP sandstone formation, α is estimated to be equal to 1 as a conservative approximation for drained deformation condition. All respective laboratory tests were conducted at S_L = 210 bar based on the calculated weight of the overburden acting on the cored sections of the wells. Since the depth difference between the deepest and shallowest cored interval is approx. 230 m, the stress difference arising from depth difference is negligible.

3.3.3. Data Processing

Petrophysical Rock Typing

Rock typing can be defined as dividing the reservoir into distinct units with characteristic petrophysical and flow characteristics [39]. Core-based Petrophysical Rock Typing methods can be classified into three separate categories:

• Methods that utilise permeability–porosity relationship and connate water saturation to some extent, excluding the so-called cut-off based methods [40];
• Methods that are based on capillary pressure data (or J-function) and measured R₃₅, e.g., Winland’s R₃₅ method, where R₃₅ is the calculated pore-throat radius at 35% mercury saturation from a mercury-injection capillary pressure test [41];
• Methods that rely on formation zone index (FZI), which is a modification of Kozeny–Carman equation, and its derivates, e.g., the spontaneous imbibition rate-driven method of FZI [42].

According to [41,42,43], the most widely used PRT methods for the classification of clastic reservoirs are Winland’s R35 method [44] and FZI-based techniques [45]. The FZI method has the advantage over the other two methods in that it allows the correlation between the micro-scale attributes and macro-scale parameters, i.e., porosity and permeability, based on the theoretical model. On the other hand, further approaches, such as Winland’s R₃₅ method, are based on empirical relationships. Since connate water saturation is unknown and no mercury intrusion porosimetry was conducted on all of the samples, we apply the FZI-based PRT technique.

According to Amaefule et al. [45], FZI-technique, Petrophysical Rock Typing (PRT) is based on grouping samples by FZI values that describe both the storage capacity (porosity) and fluid flow capacity (permeability) of the reservoir rock. This approach entails the clustering of different lithofacies of similar internal textural grain–pore compositions and petrophysical properties [46].

FZI is calculated from Reservoir Quality Index (RQI) in $\mu$m and normalized porosity (φ_Z) using the formulas below [47]:

$$\mathrm{FZI}=\frac{\mathrm{RQI}}{{\phi }_Z}$$

(2)

$$\mathrm{RQI}=0.0314\sqrt{\frac{k}{{\phi }_e}}$$

(3)

$${\phi }_Z=\frac{{\phi }_e}{1-{\phi }_e}$$

(4)

where

• ${\phi }_e$—effective fractional porosity is the ratio between pore volume and grain volume
• k—permeability in mD

Statistical Methods

In order to reduce non-normality of the dataset, permeability and grain size data were transformed to a logarithmic scale. In the case of He permeability, a base ten logarithm of the values was used. Grain size data were transformed to phi scale as proposed by Krumbein with the following formula [48]:

$$\phi =-lo{g}_2\left(d_g\right)$$

(5)

where d_g is grain diameter in mm unit.

In order to test the normality of the variables, Shapiro–Wilk tests were performed.

A Kruskal–Wallis non-parametric hypothesis test was used to investigate statistical differences between samples according to the horizontal and vertical orientation. The null hypothesis of the test is that the mean ranks of the groups are the same [49].

Petrophysical Rock Typing is based on clustering samples into groups based on their FZI values. In order to find these, mixture analysis was performed which estimates the parameters of at least two normal distributions by maximum likelihood approach in PAST software. The FZI values are divided into classes with normal distribution as a result of the non-hierarchical clustering method [50]. The optimal number of Petrophysical Rock Types was determined by Akaike Information Criterion [51]. A minimum value of AIC indicates the number of groups that produces the best fit without overfitting [50].

For statistical data analysis, IBM SPSS Statistics 29 [52] and PAST 4 data package [50] were used. For each identified PRT group, descriptive statistical parameters were calculated.

We note that the descriptive statistical properties of petrophysical parameters as well as textural petrographic properties, i.e., grain density, grain diameter, clay, and silt, as well as sand content, are discussed jointly in Section 4.2.

4. Results and Discussion

4.1. Pre-Petrophysical Thin Section Analysis

The thin section analysis reveals that the grey-light grey sandstones are characterised by well to very well-sorted grains (Figure 9). The grain size ranges from very fine to medium, but dominantly fine, and the grains are subangular to very angular, with low sphericity in morphology. It mainly consists of quartz, feldspar (K-feldspar, plagioclase), mica (muscovite, chloritized biotite), and carbonates with minor grenades and opaque minerals (coalified plant fragments, hematite), as well as zircon, apatite staurolite, and tourmaline as accessory (Figure 9). The sandstones are mainly poorly cemented by carbonates (calcite, dolomite) and clay minerals (sericite, montmorillonite, kaolinite, and illite). The micritic calcite cement occurs only in patches and narrow bands. Weak textural orientation is observed which is reflected by the presence of oriented mica plates.

The dark grey argillaceous marl and siltstone appear as massive units and as alternations of marl and siltstone laminae. In fine grained marls and siltstones, the amount of mica is significant. The darkish colour is the result of an increased amount of coaly plant fragments and clay minerals.

4.2. Petrophysical Measurement and Analysis Results

Table 3. Descriptive statistical values of 121 tested petrophysical and textural parameters.

Parameter	Unit	Mean	Median	Std. Deviation	Minimum	Maximum
Depth (RHO)	(m)	1882	1938	93.5	1741	1971
Grain density (RHO_grain)	(g/cm3)	2.710	2.709	0.015	2.683	2.750
Porosity (PHI)	(-)	0.27	0.29	0.06	0.06	0.34
Permeability (K_Klink)	(mD)	698	596	599	0.001	2157
FZI	(-)	0.315	0.332	0.176	0.003	0.615
P-wave velocity (VP)	(m/sec)	2878	2823	416	2070	4921
S-wave velocity (VS)	(m/sec)	1840	1791	284	1477	3193
Median grain diameter (d50)	(micron)	119	124	56	16	221
Clay + fine silt content (CLAY+FSILT)	(%)	16.9	14.7	8.2	9.5	49.2
Sand content (SAND)	(%)	42.3	49.7	20.9	0.0	70.1

summarises the descriptive statistical parameters of measured petrophysical parameters regarding their mean, median, standard deviation, and minimum and maximum values based on 121 samples. Since these characteristics may be biased by other factors, e.g., sampling location, rock fabric, and heterogeneity, these are classified on an unsupervised basis to reveal possible correlations. Figure 10 shows the histograms of porosity, permeability, and the calculated FZI values. These histograms exhibit multimodal distributions.

We tested the dependence of porosity and permeability on sample orientation (horizontal and vertical) using a Kruskal–Wallis non-parametric hypothesis test. The test results indicate that these parameters are independent of orientation. Therefore, Petrophysical Rock Typing was applied to the whole dataset with no respect to sample orientation.

Based on unsupervised clustering of FZI values, samples were classified into four different groups (Figure 11), where group 1 has the lowest FZI and group 4 has the largest one. These groups of samples are interpreted as Petrophysical Rock Types (PRT) that share similar reservoir characteristics.

4.3. Characterisation of Petrophysical Rock Types

Figure 12 shows the box plots of petrophysical and textural parameters for each classified PRT group. The descriptive statistics for each identified PRT are summarised in

Table 4. Descriptive statistical values of each identified petrophysical rock type (PRT). The analysed parameters are defined in Table 3.

Statistical Parameter	PRT	RHO_grain	PHI	K_Klink	FZI	VP	VS	d50	CLAY+FSILT	SAND
		[g/cc]	[%]	[mD]	[-]	[m/sec]	[m/sec]	[micron]	[%]	[%]
Mean	1	2.725	13.34	0.36	0.017	3347	2264	36	32.0	7.8
	2	2.718	25.48	61.9	0.131	2840	1867	65	18.9	22.5
	3	2.709	29.53	649	0.340	2746	1720	128	14.3	48.3
	4	2.698	29.46	1499	0.527	2862	1805	173	11.5	59.9
Std. Deviation	1	0.015	5.30	0.77	0.015	763	456	19	10.7	11.7
	2	0.013	2.43	39.9	0.044	222	77	22	3.8	13.4
	3	0.014	1.70	283	0.055	213	127	37	1.8	11.5
	4	0.009	2.42	419	0.048	258	144	33	1.7	6.6

. The figure shows that almost all tested parameters exhibit a clear dependence on the PRT group, i.e., FZI value. However, the dependence is moderate for ultrasonic wave velocities. Furthermore, higher PRT is associated with higher porosity and permeability as well as grain diameter and sand content, but with lower grain density. Regarding porosity and permeability, a clear distinction is visible between PRT 1 and PRT 2, i.e., mean porosity of 11 % versus 26%, as well as 1>> mD versus 90 mD. On the other hand, PRT 3 and 4 exhibit only a slight difference with respect to porosity. The difference between these FZI values, therefore, is associated with permeability contrast.

Regarding textural parameters, larger PRT values show an inverse relationship with grain density and clay content, as well as larger median grain size. Generally, smaller PRT values can be associated with clay and clayey appearance, while PRT 3 and PRT 4 resemble sandstone characteristics.

4.4. SEM and XRD Analysis of Petrophysical Rock Types

The results of XRD mineralogical analysis indicate that the amount of quartz does not vary between different Petrophysical Rock Types. The number of carbonate minerals (including calcite and dolomite) and phyllosilicates (montmorillonite, kaolinite, and muscovite) decreases with increasing PRT group number (Figure 13).

Scanning electron microscope analysis confirms that Petrophysical Rock Types characterised by higher FZI have a bigger average grain size (Figure 14a–c). Tangential (point) contacts of sandstone grains indicate a low level of compaction. With decreasing grain size, a higher amount of clay minerals can be observed. Authigenic clay minerals derived from weathered feldspars can reduce the initial porosity and permeability due to blocking of pore throats (Figure 14d).

4.5. Discussion

Based on the laboratory results, the petrophysical parameters, i.e., porosity and permeability, can be related to textural parameters, i.e., grain size, and clay and sand content. The textural characteristics show that primary rock textural features are not disturbed, which is a hint towards a low level of diagenetic processes, such as compaction, mineralization, and cementation. Therefore, the petrophysical properties, including Petrophysical Rock Types (PRT), can be associated with the depositional processes and textural features of the samples.

Table 4 shows that samples belonging to PRT-4 classes exhibit the highest porosity, approx. 30%, and a mean permeability of approx. 1400 mD. These sandstones can be described as clean sands with large grain size and low clay and fine silt content, approx. 10%. Sandstone samples of the PRT-3 group also show a high porosity of approx. 30%, but a lower permeability, around 650 mD. The petrographical analysis of the core samples from this class indicates smaller grain size and higher clay content. Moreover, clay content can be associated with authigenic clay minerals that are grown due to the weathering of feldspars. These clay minerals can not only reduce permeability, but they can be a potential source for fines migration eventually leading to injectivity decline. Samples assigned with PRT classes 2 and 1 show much lower mean permeability, 60 mD and 0.36 mD, as well as porosity, 25 and 13 %, and sand content compared to the previous classes. Consequently, rock samples of PRT classes 1 and 2 are more inclined to fines migration than the other classes.

These findings can be applied to field scale to make recommendations for selecting screen intervals with low potential for fines migration prior to injection or production operation. Figure 15 shows the screen intervals I to V, core sections, and the PRTs, as well as gamma ray (GR) logs of the lower sampling intervals in wells SZT-1 and SZSZT-IX (Table 2).

In well SZT-1, three out of five perforation intervals are sampled by cores. In well SZSZT-IX, core intervals are located above the perforated intervals. Regarding well SZT-1, perforation interval I of 9 m length is dominated by PRT-3 classes with few PRT-2 samples. Perforation interval II with a length of 3 m shows samples with PRT-4 class and perforation interval IV of 6 m length is associated with samples from PRT-3 and 4 classes. According to this comparison, perforation interval IV is the best candidate for sustainable reinjection. Regarding screen intervals III and V, III is less preferred due to its short zone length of 3 m, while screen interval V can be also a good candidate for reinjection operation based on the shape of the GR log.

Concerning well SZSZT-IX, it can be noted that the core interval between 1835 and 1840 m implies ideal conditions for reinjection or production operation with lower potential for injectivity or productivity problems, as the interval is associated with samples of PRT-3 classes.

Regarding processes resulting in injectivity or productivity decline in Upper Pan-nonian sandstone reservoirs, in other reservoirs with similar geological settings in Hungary, further possible mechanisms are considered as well. These may include clogging due to water–rock interaction, lack of continuous flow paths in the reservoir, and biofilm production ([6,11,20]). The investigation of these processes should be the focus of future research. Nonetheless, Szanyi et al. [3] report that in the Szeged geothermal system, in the proximity of Szentes Geothermal Field, productivity issues related to fines migration occur frequently, which may be treated by production with high flow rates.

We note that the application of the proposed methodology on GR logs for uncored intervals in both wells using machine learning techniques (e.g., [42,43]) is a subject of future research. It must be also pointed out that in our study, temperature differences between cold injection water and hot reservoir fluid during reinjection are not considered. We expect that introducing this effect may play an important role in coupled hydro-mechanical processes in wellbores drilled in unconsolidated reservoirs, e.g., fines migration due to formation damage, since fluid density and viscosity are strongly controlled by temperature [53]. The numerical study conducted by Zhang et al. 2022 [54] shows that hydraulic gradient is one of the major controlling parameters in fines production. Thus, the investigation of coupled thermal-hydraulic-mechanical analysis of near-wellbore fines migration should be the focus of future investigation.

5. Conclusions

In this study, we conducted a series of petrographical and petrophysical laboratory experiments on 121 samples of Upper Pannonian sandstone formation obtained from two exploration wells at the Szentes Geothermal Field. The goal of the study is to provide experimental data and to gain a better understanding of the formation characteristics that control injectivity and productivity issues in Upper Pannonian sandstone layers.

Based on hydro-mechanical properties of the tested rock samples, these can be classified into four representative Petrophysical Rock Types that share distinct petrophysical, textural, and hydraulic characteristics, i.e., two with higher clay content and two with higher sand content. Although sand layers are ideal for reinjection operations, one of the sandy rock types is characterised by the presence of authigenic clay that may migrate during fluid flow, resulting in injectivity decline. Consequently, the proposed methodology can be applied for identifying sand intervals with lower potential for formation damage.

The results imply that fines migration due to formation erosion is one of the key processes that must be better understood and controlled in order to mitigate injectivity issues related to the unconsolidated Upper Pannonian sandstone reservoir at Szentes Geothermal Field. However, other processes, such as mineral precipitation due to water–rock interaction processes and microbial activity, may be also considered.

Future investigation should include experimental characterisation of formation damage, including water–rock interaction tests, critical flow velocity measurements, and fines migration analysis under reservoir conditions. Furthermore, temperature effects arising from the injection of cold water into hot formation should be also studied in detail, e.g., in terms of a coupled thermal-hydraulic-mechanical numerical analysis of fines migration in wellbores in unconsolidated reservoirs.