Structural and Parametric Models of the Zalezcze and Zuchlów Gas Field Region, Fore-Sudetic Monocline, Poland – An Example of a General Static Modeling Workflow in Mature Petroleum Areas for CCS, EGR or EOR Purposes

— Za łę cze and Ż uchlów are strongly depleted natural gas ﬁ elds in aeolian sandstones of the Rotliegend, located in the central part of the Fore-Sudetic Monocline. A set of three static 3D models was generated to check the possibility of CO 2 injection for Enhanced Gas Recovery (EGR) and to check the safety of storage by means of geomechanical modeling: one regional model (ZZA) and two local models – the ﬁ rst for Za łę cze (ZA) gas ﬁ eld and the second for Ż uchlów (ZU) gas ﬁ eld. The regional model is composed of 12 stratigraphic complexes (zones) from the base of the Rotliegend to the ground surface. The local models comprise only the three lowermost complexes: ﬂ uvial deposits of the Rotliegend, aeolian sandstones of composé de 12 complexes stratigraphiques (zones) depuis la base du Rotliegend jusqu ’ à la surface du sol. Les modèles locaux comprennent uniquement les trois complexes les plus bas : dépôts ﬂ uviaux du Rotliegend, sables éoliens du Rotliegend (Réservoir I) et calcaire basal de Zechstein, Ca1. Les éléments clés de la procédure de modélisation comprennent : le contrôle qualité (CQ) des données, l ’ interprétation des paramètres manquants nécessaires pour la modélisation statique et leur intégration dans un géomodèle. L ’ approche a été élaborée pour produire des modèles régionaux et locaux convergents. Le modèle statique régional représente un cadre pour un modèle géomécanique régional. Les modèles locaux permettent des simulations dynamiques et une modélisation géomécanique locale. La méthodologie présentée pourrait être utilisée moyennant quelques adaptations à la géomodélisation de nombreux gisements de gaz et de pétrole matures.


INTRODUCTION
This paper presents an outline of the methodology and workflow applied for the creation of three complementary 3D static geomodels, for testing the possibility of CO 2 injection for underground storage (Carbon Capture and Storage, CCS) and Enhanced Gas Recovery (EGR), as well as for the security assessment of these processes by means of geomechanical modeling. The models cover an area surrounding two big depleted gas fields, Załęcze and Żuchlów, and many minor gas accumulations.
Załęcze and Żuchlów gas fields are located in South-Western Poland, on the Fore-Sudetic Monocline (FSM). The aeolian and fluvial sandstones of the Upper Rotliegend are the major reservoirs in this area. They are covered with Zechstein Basal Limestones (Ca1), which is a reservoir in the region. Zechstein deposits above the Ca1 form the caprock and an effective seal of the Rotliegend reservoir. The anhydrites and salts of the PZ1 Cyclotheme form the direct seal ( Fig. 1, 2).
Exploration in the Załęcze and Żuchlów area allowed the discovery of several gas accumulations. Załęcze and Żuchlów were the biggest of them (Fig. 3). Exploration and development work was finished before 1989, with 231 completed wells. They had relatively good coring of the reservoir section and very abundant but generally poor geophysical logging. The region is not covered by 3D seismic surveys and only partly covered by one 2D survey from the seventies to nineties. In this situation, the crucial effort was to combine and optimize the input data both for structural and property modeling, and then to elaborate modeling workflows, allowing obtaining concordant regional and local static models. Generally staying in the mainstream methodology, modeling was completed with original tailor-made workflows, fitted to the multipurpose use of the model and the quantity and quality of the input data. The basic regional model (ZZA), built for the whole research area, is composed of 12 units (zones) extending between the base of the Upper Rotliegend and ground surface. The number and type of layering of zones are linked to their genesis, thickness, internal variability and further processing destination. The entire structural regional model is composed of about 11 9 10 6 grid cells. The 3D grids of the Rotliegend zones comprise models of facies, superimposing local facies variability on the general sedimentary pattern reflected by zonation. Property models were based mainly on log data. They were created for further use in dynamic modeling and to illustrate the variability of the basic parameters governing gas field quality: effective porosity (log-based and core-based); shale volume (V sh ); permeability (log-based and core-based); relative permeability (K w , K r ); and water saturation (S w ). For further geomechanical modeling preliminary log-based parametric models were elaborated, including acoustic impedance, density RHOB and velocity V p .
The ZZA is the framework of the two local models covering only the three lowermost complexes. The first of the local models covers the big Załęcze and marginal Wiewierz gas field (ZA) area and the second corresponds to the Żuchlów (ZU) gas field area. The local models were extracted from the regional one, testing two methodologies in order to get the best possible concordance with the flow simulator during dynamic modeling of CO 2 injection for EGR (Klimkowski et al., 2015).
The modeling workflow presented in this paper can be adopted for EOR, EGR or CCS models in many depleted oil/gas fields of Poland, and more generally of Eastern Europe and Russia, where the methodology of exploration and development was in general very similar before 1990. Working on these kinds of mature fields, geomodelers face the problem of unequivocal and effective integration of very numerous and low-quality data. Lack of high-quality data focuses the workflow on more traditional modeling methods, demanding first of all geological expertise.

GEOLOGICAL SETTINGS
The Załęcze -Żuchlów area is located in the southern part of the Fore-Sudetic Monocline ( Fig. 1) (Narkiewicz and Dadlez, 2008;Karnkowski, 2008;Żelaźniewicz et al., 2011), which is a tectonic unit comprising the Permo-Mesozoic and covers the Paleozoic Platform (PP). The PP in this area contains the pre-Carboniferous basement, overlain by the Carboniferous fold-and-thrust Variscan belt (Żelaźniewicz et al., 2011).
Sedimentation of the Permo-Mesozoic rocks covering the FSM area was in general continuous from the Rotliegend up to the Upper Cretaceous (Marek and Pajchlowa, 1997). The interval region was situated at that time on the southern flank of the Polish Basin, known as the Mid-Polish Trough (MPT), which in turn was part of the Permian-Mesozoic system of epicontinental basins of Western and Central Europe, overlapping to a great extent the area of the Southern Permian Basin (SPB). The Polish Basin was the eastern part of the SPB (Ziegler, 1990;Kiersnowski et al., 1995;Krzywiec, 2006;Jarosiński et al., 2009;Gast et al., 2010;Pharaoh et al., 2010). Its development resulted from longterm thermal subsidence, which started in the Permian and lasted until the Late Cretaceous. It comprised three major pulses of extension-related accelerated tectonic subsidence during Late Permian to Early Triassic times, in the  Oxfordian to Kimmeridgian, and in the Early Cenomanian (Dadlez et al., 1995;Stephenson et al., 2003;Krzywiec, 2006;Sowiżdżał et al., 2013). Partial inversion, which took place in the FSM area during the Late Cretaceous under a NE-SW directed transpressional regime (Deczkowski and Gajewska, 1980;Lamarche et al., 2002;Jarosiński et al., 2009) resulted in complete erosion of Cretaceous, partial to complete erosion of Upper, Middle, Lower Jurassic and partial erosion of Triassic strata (Fig. 2) (Górecki, 2006a, b). Cenozoic strata deposited north of the Sudetes are essentially undeformed, excluding active salt-induced structures and narrow grabens, 1-2 km wide, several hundred meters deep and several kilometers long (Jarosiński et al., 2009).
Tectonics of the ZZA comprise the pre-Permian fault systems trending NW-SE (e.g. the Middle Odra and Dolsk Fault Zones), SW-NE and W-E. The W-E system was reactivated during the Triassic -Early Jurassic tensional regime, resulting in the development of negative flower structures within the Mesozoic cover (Deczkowski and Gajewska, 1980;Grocholski 1991;Kwolek, 2000Kwolek, , 2003Jarosiński et al., 2009). These structures, although rooted at the base of the Zechstein, create separate Mesozoic fault systems (Markiewicz, 2007). The youngest fault system is closely related to Mesozoic tectonics. This Cenozoic graben system developed during the latest Eocene and Oligocene times. The grabens propagated northwards into the Polish Lowlands at the end of the Oligocene and remained active until Mid-Miocene times (Jarosiński et al., 2009;Pharaoh et al., 2010).
The ZA and ZU gas fields belong to the Carboniferous Total Petroleum System (CTPS) (Magoon and Schmoker, 2000), one of the most important European groups of petroleum provinces (Glennie, 1998;Karnkowski, 2007;Pletsch et al., 2010;Botor et al., 2013) which, in the area of the FSM, is a marginal part of the Southern Permian Basin. The aeolian and fluvial sandstones of the Upper Rotliegend form the major reservoirs in this area.
Deposits of the Rotliegend group, represented in the Polish part of the SPB by volcanic and terrigenic series, are organized into the number of formal lithostratigraphic or allostratigraphic subdivisions (Pokorski, 1981a(Pokorski, , 1988Karnkowski, 1994Karnkowski, , 1999bKiersnowski, 1997Kiersnowski, , 1998Kiersnowski and Buniak, 2006). An attempt was made recently to unify these subdivisions and correlate them with those in the German part of the SPB (Gast et al., 2010). Additionally, simple subdivision of the Rotliegend into non-reservoir Autunian and reservoir Saxonian, accompanied by paleogeographic division, is commonly accepted in prospecting practice. The number of paleogeographic units was distinguished in the Polish part of the Rotliegend basin. The Central and Silesian basins form the largest two units (Pokorski, 1981b(Pokorski, , 1997Karnkowski, 1994;Kiersnowski, 1997;Kiersnowski and Buniak, 2006;Papiernik et al., 2010aPapiernik et al., , 2012. The ZZA is located close to the center of the latter. Rotliegend deposits in the Silesian Basin are developed in 50-60-km-wide complexes. This zone is bordered from the NE by the Brandenburg-Wolsztyn-Pogorzela High and from the SW by the Fore-Sudetic block (Fig. 3). The basin was created due to tectonic movements during Late Carboniferous and Permian times, strongly influencing deposit thickness and facies distribution. The paleo-relief was flat during the Upper Rotliegend time. This caused facies unification of overlapped Upper Rotliegend deposits. The Upper Rotliegend deposit thickness is generally equal within the Silesian Basin, but locally exceeds more than 400 m in places of most subsiding tectonic blocks due to still active tectonic movements. These deposits are represented mostly by aeolian sands occupying the central part of the Rotliegend Silesian Basin (called Southern Erg) and by coarser fluvialalluvial deposits in the basin margins and the deeper part of the profile (Fig. 3). Several large dune fields originated in the central part of the Rotliegend Silesian Basin in the last stage of deposition. These dune fields were partly preserved during transgression of the Zechstein Sea and, as a result, formed geomorphological traps for gas accumulation, housing at least 20 gas fields. Załęcze and Żuchlow belong to the largest among them (Kiersnowski and Tomaszczyk, 2010;Papiernik et al., 2012).
Rotliegend reservoirs are covered with Zechstein Limestones (Ca1), which in the investigated region essentially are not producing reservoirs, although they produce gas south of the ZZA. In the Brandenburg-Wolsztyn -Pogorzela Height, the Basal Limestone (Ca1) is developed in reservoir facies. Zechstein deposits form the caprock and the effective seal of the Rotliegend reservoir above the Ca1. The anhydrites and salts of the first Zechstein Cyclotheme, PZ1, form the direct seal of the CTPS. They are superimposed by the main dolomite carbonates of the PZ2 Cyclotheme, gas-and oil-bearing deposits being the source rock and the reservoir of the Zechstein Petroleum System (Karnkowski, 1999a;Pletsch et al., 2010). Within the area of the Fore-Sudetic Monocline, Ca2 can be treated as the secondary reservoir for CCS purposes. Its caprock is built of anhydrites, salts and shales of PZ2, PZ3, and partly PZ4 Cyclothemes. In the investigated area, PZ2-PZ4 Cyclothemes are reduced and they do not contain all units marked in Table 1 (Peryt et al., 2010).
The stratigraphic column above the Zechstein comprises clastic deposits of the Lower and Upper Bunsandstein (Tp1+Tp2), carbonates of the Middle Triassic (Roet -Tp3 and Muschelkalk -Tm) in the ZZA.
Erosional remnants of the Upper Triassic are present only close to the Northern margin of the ZZA. The Tertiary and Quaternary deposits form the topmost part of the ZZA profile. Their thickness is up to a few tenths of meters.  Figure 2 Regional section through the Fore-Sudetic Monocline -ca. 15 km SE from Załęcze gas field without Cenozoic deposits (based on Górecki, 2006a, Doornenbal andStevenson, 2010). Geological reservoir modeling workflows are now commonly used in the oil industry, from the construction of complex structural grids to their filling with petrophysical properties, honoring well data. Geostatistical simulation methods are classically used at that step as they produce realistic images of the geological complexity and internal distribution of reservoir heterogeneities, honoring the data. Geostatistical approaches also allow the quantification of the uncertainty associated with the results obtained. As also already noted by Journel and Alabert (1990), and illustrated in their paper, data integration provides a major contribution to reservoir modeling geostatistics. Moreover, this approach allows the integration of secondary data, such as qualitative conceptual geological models (Ravenne et al., 2002) or quantitative seismic data (Lerat et al., 2006;Doligez et al., 2007).
In this study, the modeling was performed starting from typical workflows for the reservoir's units and its overburden (Doligez et al., 1999(Doligez et al., , 2007Cosentino, 2001;Dubrule, 2003;Labourdette et al., 2005;Zakrevsky, 2011). However, according to the available data, the stratigraphic framework was made the focus. The stratigraphic division was defined in order to allow the optimum representation of reservoir variability and geomechanical parameter uniformity in the reservoir as well as in the overburden.
Firstly, a detailed regional model of the reservoir and its overburden was elaborated for the whole ZZA area. Due to the different quality and quantity of data, different workflows of property modeling were launched for the main reservoirs and their overburden.
The horizontal and vertical variabilities of the reservoirs were estimated through laboratory data (effective porosity measurements -PORO, and total permeability -PERM) as well as log data comprising shale volume (V sh ), porosity (u), total permeability (K ), water permeability (K w ), gas permeability (K g ) and water saturation (S w ). We modeled all listed properties using 10 stochastic realizations of Gaussian Random Function simulations. Finally, base-case models of properties were estimated as an arithmetical or geometrical average of these 10 realizations depending on the distribution type of the parameter.
The models calculated with this workflow used hybrid (deterministic and stochastic) solutions. They offer smoother results than individual stochastic realizations; however, they reflect the uncertainty of results away from data better than the pure deterministic solutions. This procedure was applied only to the petrophysical data within the reservoir zones, where the abundance of data, their spatial distribution and relatively higher confidence allowed higher control of the stochastic models by means of (local) shorter-range variograms, variable distribution and conditioning data for the resulting models.
The property modeling workflow was performed to estimate the facies distribution in the Upper Rotliegend and Zechstein Limestone (Ca1) complexes, as well as the controlling parameters of the overburden such as V sh and u. The process comprised data analysis, variogram construction and modeling with the use of indicator kriging for facies and simple kriging, reproducing relatively simplified but reliable trends of variability for V sh and u.
The study area is relatively large and the spatial well distribution is very uneven, with clusters of wells within gas field areas. Such a distribution together with the lack of geo-steering data, such as e.g. AI or depth 3D seismic cubes, implied inevitably huge marginal errors among the regional and local models calculated for the ZZA separately. To avoid these discrepancies, it was decided to generate a regional model that was as detailed as possible and then to extract the local models from it. The reservoir models of the ZU and ZA gas fields are thus parts of the regional model, comprising its three lowermost zones. Two methods of their extraction were tested. First, successful tests of populating the local simple grids ZA and ZU with the variables from the regional grid using a process of structural and parametric upscaling were performed. The quality of these models was satisfying, with some minor errors of parameter assignment to the zones. However, finally these models were rejected, mostly due to practical problems with proper import to external dynamic simulation software.
The second, more convenient method of extraction is the procedure of simple copying of selected zones and cells of the regional model. It allowed us to maintain full concordance of the geometry, and parametric model variability. This simple procedure is especially useful, as it does not introduce any geometrical errors, which could impede the operation of the dynamic simulation software.

Input Data
The geometrical framework of both the local and regional models was based on structural maps constructed as 2D grids. These maps were only partly based on interpretation of the relatively new 2D seismic data, reprocessed in 2001. At least half of the area of interest, especially the Żuchlów zone, is not controlled by digital seismic data of good quality. To fix the problem of lacking control, the regional 2D grids representing structural and thickness variability were used and corrected locally with abundant well stratigraphic information, and horizontally refined. To solve problems of partial coverage of the AOI with seismic data, we created a set of complex maps leading to a patchwork of local seismic maps and regional models (Fig. 4) as well as a series of intermediate surfaces generated as a result of map operations, refinement and well top corrections.
The most detailed cartographic input data comprised seismic maps (2D grids) of four seismic horizons: the Zechstein base, the Z2 reflector (the top of the basal anhydrite A2), the top of the Zechstein and the top of the Upper Bundsandstein, Tp2. They were fitted to the well tops of the Saxonian, Main Dolomite (Ca2), Zechstein's red pelite (IP) and Upper Bundsandstein (Tp2). These materials were too sparse to create a framework ensuring the required resolution.
To increase the stratigraphic resolution of the model, regional trend maps in the form of 2D grids were used. These maps included: a structural map of the Zechstein base (Saxonian top) (Fig. 4), a thickness map of the Zechstein Limestone -Ca1 , a thickness map of the Upper Rotliegend (Saxonian) , a thickness map of the Upper Rotliegend fluvial facies deposits , a thickness map of the Upper Rotliegend aeolian facies deposits , a structural map of the Zechstein top, resampled high-resolution topography 2D grids (Tomaszczyk, 2007) as the topmost surface of regional structural models. Stratigraphic depth information from 231 wells, located within the AOI, was used either to correct seismic/regional 2D grids locally or to generate lacking surface (Ca2 base) during structural 3D modeling (Fig. 5).

Fault Model
The fault model comprises a very general interpretation of fault geometry based on 2D seismic and regional maps. The spatial range of the original fault model reaches the area outside the regional static model (Fig. 6). Faults in the Zechstein substratum run mostly in the SW-NE direction. Some smaller faults developed in a roughly perpendicular direction (SW-NE) were observed at the base of the Upper Rotliegend. Faults developed within the Rotliegend deposits do not cause compartmentalization of gas fields, as they have very Area of regional 3D model Local seismic map Figure 4 Location of the regional 3D model and input seismic map on the background of a regional 2D grid of the Zechstein base The faults observed above the Zechstein form the SW-NE grabens and display throws up to 200 m (Fig. 6, 7). They were generated in a transpressional regime as a result of reactivation of the Paleozoic fault system. However, the two systems' continuity is disrupted by the plastic salts and anhydrites of the Zechstein (e.g. Załęcze gas field). The geometry of the grabens is very poorly controlled but correlation of seismic and well data shows that in some places, e.g. above Załęcze gas field, local-scale diapirism can generate some local inversion which is not delineated by 2D seismic data.
The precision of the fault model is not very high as a consequence of the low quality of the input data. Another problem concerning fault modeling is the occurrence of the two independent fault systems.

Structural Model
The most complete structural regional model was adapted to geomechanical modeling requirements. It contains 12 horizons whose estimations are based on structural and/or thickness maps and corrected on suitable well top positions: topography,

Model Zonation and Layering
Using the listed horizons, the regional model was divided into 12 complexes or zones according to their geological genesis, facies or functionality for modeling (Fig. 7). In the Zechstein part, zones are the equivalents of lithostratigraphic units presenting uniform lithologies. The division of the Zechstein is detailed from the base of the Ca1 up to the top of the Ca2. The PZ2, PZ3 and PZ4 cyclothemes were divided into two zones. More detailed modeling was not possible due to the low thickness of lithostratigraphic units and the poor quality of well picking. The zonation in the reservoir part follows the general vertical lithological variability.
The three lowermost zones comprise the reservoir section of the geological profile for further simulation purposes. Identification of these zones was based on regional sedimentology and rather detailed core and log facies analyses. Within the crucial reservoir zone, the Saxonian aeolian, an optimum layer thickness was established as 10 m, based on the vertical variability of the u logs. A layering following the base was applied to reflect the internal architecture of the mega-dunes in the ZZA. The same layering type was also used for the Ca1 carbonate zone, where the layer thickness was set to four meters to account for petrophysical variability.

B. Papiernik et al. / Structural and Parametric Models of the Załęcze and _
Zuchlów Gas Field Region, Fore-Sudetic Monocline, Poland -An Example of a General Static Modeling Workflow in Mature Petroleum Areas for CCS, EGR or EOR Purposes for geomechanical modeling. Each zone includes similar terrigenic lithological successions. Layering of non-reservoir complexes is not very dense. It was determined with the trial-and-error method to optimize the model size and to avoid over-averaging of petrophysical parameters, compared with log scale variability. For all these zones, proportional layering starting from the base with minimum thickness was applied. Table 2 summarizes the zonation and layering of the regional model.
The local reservoir models of the ZA and ZU fields have the following geometrical statistics: the ZA model comprises 114 9 118 9 72 nodes, 72 geologic layers and 21 faults within 3 zones; the ZU model comprises 91 9 76 9 72 nodes, 72 geologic layers, 5 faults, and 3 zones.

PROPERTY MODEL
Property and facies modeling was performed within the regional structural framework and then the results were transferred into the local reservoir ZA and ZU models.
Modeling of the overburden of the Rotliegend was very simplified, due to the lower resolution of geomechanical models, and the amount of available data (only logs) in this interval.

Input Data for Facies Modeling
In the investigated region, facies logs from joint geophysical logs and core descriptions were relatively sparse and available for only six wells in the field area and four wells close to the grid boundary (Fig. 8). Despite the small quantity of logs, the facies model of the Saxonian aeolian reservoirs seems to be highly reliable as the variability of sandstones  in this complex is rather low. On the other hand, the lowermost fluvial zone of the Saxonian is to some extent conceptual due to the smaller number of data. Regional facies interpretation was based on complex core and log interpretation  and was supplemented with some additional GR log interpretation within the ZZA. Both interpretations resulted in discrete lithofacies logs comprising aeolian (1), fluvial (2) and marginal playa (3) facies. The latter occurs only at the outskirts of the ZZA and is negligible.

Input Data for Property Modeling
Input data consisted only of laboratory data and log interpretations and comprised: laboratory measurements of porosity and permeability (191 wells); geophysical logs of shale volume (V sh ), porosity (u), total permeability (K ), water saturation (S w ), permeability for the water phase (K w ), and permeability for the gas phase of the Upper Rotliegend and Zechstein Limestone, Ca1 (K rg ) (45 wells); logs of shale volume (V sh ), porosity (u), water saturation (S w ), sand volume (V sand ), halite volume (V halite ), anhydrite volume (V anhydrite ), limestone volume (V limestone ) and dolomite volume (V dolomite ) as well as Velocity (V ), density (RHOB), and Acoustic Impedance (AI) logs of the reservoir and its overburden (26 wells located at the ZZA and close to its margins).

Facies and Property Model
The result of facies and property modeling is basic parameter models. Their spatial variability was approximated in the regional ZZA model, using a mixed methodology of estimation, individually fitted to data variability in the individual zones of the model. Table 3 summarizes the methods and parameters used to populate the regional model with these properties. In addition, some preliminary properties suitable for the next phases of modeling were estimated, including bulk density (RHOB) and Acoustic Impedance (AI). These models were later improved and/or replaced with more sophisticated ones elaborated by geomechanical modelers.

Facies Models
Facies models were created for reservoir zones, as only there were suitable data available. In general, the stratigraphic zonation of the overburden was considered as lithologically uniform stratigraphic complexes.
Within the reservoirs, the facies model is also simple as a consequence of the genetic zonation and rather low resolution of input data. Four basic facies systems were defined: aeolian sandstones (code 0); fluvial and alluvial deposits (code 1); marginal playa deposits (code 2); and undivided carbonates (code 3) (Fig. 9).
The procedure of facies modeling of the three potential reservoir zones is summarized in Table 4. Carbonates of the topmost Ca1 -Carbonates (Reservoir I) zone in the ZZA area represent open basin deposits displaying low thickness and very poor reservoir properties (Peryt et al., 2010). As a result of such uniformity, the whole zone was assigned as carbonate facies (code 3).
The sandstones of the Upper Rotliegend building the Saxonian aeolian zone (Reservoir I) are productive rocks within the ZZA. According to data analysis, this complex is dominated by aeolian sandstone at 95.5% (code 0). The intercalations of fluvial deposits (code 1) were recorded within data only in its basal parts, close to the ZZA margins (Tab. 4). Facies in this zone were modeled using indicator kriging with isotropic exponential variograms with a range of 5 000 m. These values were adopted arbitrarily as data analysis does not work properly with such a homogeneous input data set.
The Saxonian fluvial zone remains in hydraulic contact with the aeolian reservoir but is very distant from the gas/water contacts of the fields existing in the ZZA. So, it can be treated as a substratum of the effective reservoir complex. For this reason, its layers display relatively high thickness. The investigated complex is composed mostly of fluvial deposits (63%) with a large admixture of aeolian sandstones (~33%). Marginal playa deposits (code 2) observed in input data (3.5%) (Tab. 4, Fig. 10) are located outside the ZZA and in practice they are negligible in the modeled area. The poor data control can be a considerable problem influencing modeling quality, as most of the wells used for facies modeling in the ZZA do not reach the base of the Saxonian. As a result of data analysis, modeling was achieved using spherical variograms strongly anisotropic for aeolian intercalations and practically isotropic for fluvial deposits. The variogram for playa deposits is in practice meaningless. The applied variograms have large ranges, as they were designed to disclose regional features of facies variability, matching sparse data and paleoenvironmental variability.

Settings of Petrophysical Modeling
Modeling of the overburden zone properties was based on rather sparse but evenly distributed data, spread over the modeled area and close to its margins. This, together with dense and precise zonation, allowed the data analysis to show uni-modal data distribution. For these sets of properties long-range distance variography was applied to disclose regional constituent variability.
According to the data number, quantity and quality, the overburden models were designated as the best case models represented by the kriging model. Figure 11 displays an example of data analysis and the variogram for RHOB in the Lower Triassic Tp1+Tp2 zone.
The well spatial definition of salt-or anhydrite-dominated zones of the Zechstein allowed us to assume a uniform distribution of properties and arbitrary (high) range of variograms as the best matching data variability for kriging.
The properties in the reservoir units, for which the number of data allowed more detailed data analyses and parameters, were simulated with Sequential Gaussian Simulation (SGS) in order to reflect local variability. Figure 12 gives an example of representative settings of parameters used in modeling of reservoir zones, displaying the data distribution and related variogram for the major axis direction.

Regional Model Examples
A direct facies control was applied to estimate the shale content distribution (V sh ) in reservoir zones. This 3D volume of V sh was in turn used in collocated co-kriging to improve the quality of porosity distribution (u, Fig. 13) 3 Reservoir zones Hybrid modelregional background from simple kriging, and very small range variograms; local models on gas field areas average of 10 SGS realizations estimated using collocated co-kriging with the geometrical parameter height above the contact 69 Height above the contact 3 Reservoir zones Geometrical modelingwithin small regions covering gas field areafeature related to GWC and the top of the reservoir 4 gas fields The water saturation (S w ) model is an example of a hybrid model incorporating stochastic and deterministic solutions (Fig. 14). Its background displaying high S w was estimated from low-range kriging, but regions covering gas field areas were calculated separately as an average of 10 SGS realizations, co-kriged with the height over the contact, and integrated into the background model.  4 000 6 000 8 000 10 000 12 000 14 000 16 000 18 000 20 000 22 000 24 000 26 000 28 000 30 000 32 000 Facies Figure 9 Facies model of the Upper Rotliegend displayed in the W-E section through the ZZA. An enormous quantity of input data and output results does not allow for a brief summary. Table 5 summarizes the set of averaged modeled parameters collected for the Saxonian Aeolian (Reservoir I) as an illustration of petrophysical parameters.

Local Models
In order to keep a full concordance between the regional and local static models, selected zones, columns and rows were extracted by a simple copying procedure to the local models. Some minor discrepancies between the regional and local scales can arise when an upscaling technique is applied, and local updates were made. The potentially weak point of this workflow is the need for using a rather powerful CPU to generate a high-resolution regional model to reflect local variability accurately.
Regardless of the extraction method, the model can be locally updated in the presented workflow.   Shale volume (V sh ) model in the ZA area.
The average variability of the reservoir parameters within the most important producing Ps-aeolian zone are displayed for the ZA area in Table 6 and for the ZU gas field in Table 7.
The two local models were extracted and transferred to geomechanical and reservoir engineering teams. They were prepared with the use of the first of the extraction methods described. Testing revealed that it ascertains better compatibility with modeling software of other brands than Schlumberger.

Figure 16
Facies model in the ZU area.

CONCLUSIONS
The work presented allows us to draw a set of conclusions and recapitulations related to the specific results of modeling in the Załęcze-Żuchlów Area and some general conclusions on the workflow presented. Regarding the specific results for the ZZA, the following aspects of the work can be highlighted: the models presented in this paper are based on a very extensive input database covering: very rich stratigraphic identification in wells (up to 231 wells); seismic identification based on low-quality 2D seismic data, and covering only part of the project area. There are no 3D surveys and no seismic inversion results in the ZZA; good well log characterization of the reservoir section of the area (up to 71 wells) and considerably worse log control of the overburden (up to 26 wells). The quality of interpretation is in each case limited, due to the quality of the logs acquired during the years 1970 to 1990; good but uneven control of core measured porosity and permeability within the Ca1 and the topmost 50 m of the Rotliegend; relatively general quantitative facies identification based on core and log data; good regional trend map coverage; the variability of the data involved a first basic step to integrate all these data into a common digital format, allowing the unification, quality check and finally improvement of the models. A great quantity of data does not correspond to their quality. Especially, lack of 3D seismic and poor coverage with 2D surveys dramatically increase the structural and tectonic uncertainty of the ZZA models. The considerable age of the geophysical log data and shortage of modern petrophysical and geomechanical core data considerably confine the methodology of 3D modeling. No sophisticated analyses based on seismic inversion and advanced log interpretation could be applied; in spite of the considerable problems, the models elaborated allowed the successful restoration of geological realism, on the regional and local scales. Following this work, the potential CO 2 storage capacity has been estimated for the Załęcze and Żuchlów sites. The possible trapping mechanisms in these reservoir sites were also evaluated. Geomechanical modeling studies were also carried out to evaluate the geomechanical effects associated with gas extraction and potential future CO 2 injection into the two depleted gas fields Załęcze and Żuchlów.
Regarding the general features of the workflow presented, the following aspects may be discussed: the key component of the procedure was to elaborate workflows of modeling allowing us to preserve concordance between the regional and local reservoir models and a regional model suitable for further modeling; the high stratigraphic resolution of the structural framework plays a crucial role in the proposed workflow. Well-defined stratigraphic identification guided by geological expertise can improve the quality of property modeling and leads to a better spatial definition of the input that should be modeled as individual complexes (zones).
To some extent, in areas of poor data quality, well-defined zones are used for data subset extraction as simple facies models. In other cases, well-defined zonation improves the model identification; data analysis and model populating schemes were designed in a way allowing low time consumption and unambiguous delineation of parameter variability in the form of the best guess model for overburden zones; the proposed methodology of modeling is matched to the expected resolution of individual zones as well as to data quantity and quality. Generally, for the reservoir zones where data control is good, the procedures of data analysis and modeling aimed to express both regional and local variability. Modeling was performed in the form of a smoothed average of 10 stochastic property models for each parameter, being a kind of hybrid between stochastic and deterministic solutions. The resulting models display some variability in the property distribution better than typical deterministic models and they are not as "noisy" as individual stochastic realizations, but the procedure is more time-and computer memory-consuming; concerning the overburden zones where input data control is bad or moderate and where an extremely high resolution of the model is not necessary, kriging-based models were used and display a smooth regional variability of petrophysical parameters. The workflows presented can be used as a starting point for geomodeling of many mature petroleum areas such as the ZZA. This example seems to represent many hydrocarbon producing areas in Poland quite well, as well as in the former Soviet Block countries, where the methodology of exploration and development applied before 1990 was very similar. Certainly, in specific sites the problems encountered will be slightly different, but in each case the geological expertise will be crucial for data integration and optimum stratification of the model.