www.geologia-croatica.hr Journal of the Croatian Geological Survey and the Croatian Geological Society 2025 | 78/2 | 151–172 | 4 Figs. | 8 Tabs. | 1. INTRODUCTION The Geochemical Atlas of Croatia (GAC) provides insights into regional and local topsoil geochemistry and indicates the importance of further research (HALAMIĆ & MIKO, 2009). The collected data of the GAC was stored in the Geographic Information Software (GIS) database and therefore offers a possi bility of comparison with other databases that contain the same type of data. Therefore, the various influences on the geochemistry of the topsoil could be analysed. The regional variability of the natural geochemical background is also bet­ ter known, allowing better understanding and identification of anthropogenic influences. The total number of sampling points for GAC in topsoil across Croatia was 2,521, sampled in a re­ gu lar grid of 5 x 5 km. Based on the results obtained, individual maps for 27 chemical elements were generated. All elements showed differences in the spatial distribution of two geo dyna­ mic units in Croatia to some extent: Dinaric­coastal (DIN) and Pannonian Basin (PAN). HASAN et al. (2020) elaborated on the geochemical signature by geology, soil, land use and land cover and region in DIN. It is sufficient to highlight here that 1,254 topsoil samples were collected to explore the factors behind the characteristic geochemical signature of the Pannonian Basin. This was done using different environmental and geological criteria that were, by default, independent of the soil geochemistry and Geochemical Map of Croatia – Key findings and legacy of geochemical mapping in the Pannonian Basin Ajka Pjanić1, Josip Halamić1, Lidija Galović1*, Zoran Peh1 1 Croatian Geological Survey, Sachsova 2, 10 000 Zagreb, Croatia; (*corresponding author: lgalovic@hgi-cgs.hr; apjanic@hgi-cgs.hr; jhalamic@hgi-cgs.hr; zpeh@hgi-cgs.hr) doi: 10.4154/gc.2025.08 Abstract The Geochemical Atlas of Croatia revealed differences in the spatial distribution of geo- chemical elements between the Croatian Dinaric coastal region and the Pannonian Basin. In the Pannonian Basin, 1,254 topsoil samples (0–25 cm) were analysed for 27 elements and pH values. This geodynamic unit is characterised by the variability of geological, pe- dological and land-use-related characteristics. To analyse the geochemical properties of soil above different lithological units, soil types and land use and land cover classes, a dis- criminant function analysis was performed using the compositional data. The study shows significant geochemical and ecological patterns, with the first two discriminant functions explaining a major proportion of the variability: in the GEOLOGY model, DF1 explains 33.31% and DF2 17.97%; in the SOIL model, DF1 explains 51.59% and DF2 20.22%; in the CLC model DF1 explains 59.85% and DF2 30.30%. The geological model distinguishes between the Quaternary sediments and the older lithological units and highlights alkaline and acidic soil conditions. The soil model shows the effects of fluvial and alluvial deposits, agricultural practices and the underlying geology on soil composition and emphasises the enrichment of essential nutrients and heavy metals in the soils. The land use and land cov- er model illustrates the anthropogenic influence on agricultural soils and the susceptibility of wetlands to heavy metal accumulation. The results obtained illustrate the complex inter- action between geology, the topsoil layer, as well as land use and land cover, providing awareness for environmental management and monitoring, and the need for further geo- chemical studies of the soils in the Pannonian Basin. involved such divisions as surface geology at the sampling point, various soil types, and descriptions of the land use. The latter was borrowed from the land cover classes described in the Corine Land Cover (CLC; CORINE LAND COVER, 2018). In contrast to the earlier investigated Dinaric region, no REGION model was created for PAN due to the homogeneity of the sub­regions within the Pannonian Basin. Thus, three divisions – GEOLOGY, SOIL and CLC – were established, providing the most effective ways of an a priori arrangement of the soil samples into many comprehensible and all­ embracing statistical groups. In the final analysis, discriminant function analysis (DFA) was employed as a method of data reduction and structuring, generating the respective models based on geochemical partitioning between the established groups. Being mathematical (statistical) by their nature, these models created appropriate structural patterns that helped explain the behaviour of the observed geochemical data in the process­form terms (STRAHLER, 1980), most thoroughly via the geochemical maps generated in GIS. The Pannonian Basin is characterised by a continental, moderately warm and humid climate (FILIPČIĆ, 2023). The older rocks in the basin are more strongly eroded, and we do not find many of their outcrops on the surface. The basin consists mainly of Quaternary formations, including loess, Article history: Manuscript received: July 12, 2024 Revised manuscript accepted: March 5, 2025 Available online: June 30, 2025 Keywords: Geochemical mapping, geology, FAO, Corine, discriminant function analysis, compositional data, Pannonian Basin, Croatia mailto:lgalovic@hgi-cgs.hr mailto:apjanic@hgi-cgs.hr mailto:jhalamic@hgi-cgs.hr mailto:zpeh@hgi-cgs.hr G eo lo gi a C ro at ic a 152 Geologia Croatica 78/2 aeolian sands, alluvium, lake and swamp sediments, which overlie the older rock formations (HGI­CGS, 2009). The pedological composition of the Pannonian Basin in Croatia is heterogeneous, which is due, among other factors, to considerable differences in subsurface lithology, climate, agricultural production, land use and water saturation (BOGUNOVIĆ et al., 1998; BAŠIĆ et al., 2007; BAŠIĆ, 2013). The Corine Land Cover dataset provides detailed infor­ mation on land use and land cover (CORINE LAND COVER, 2018). This dataset offers insight into various land cover features such as forests, grasslands, wetlands, urban areas and water bodies. The CLC is used in spatial data studies for the management of environmental and urban areas and natural resources, among others. The GAC dataset, which includes 28 variables of topsoil, was analysed using the multivariate statistical method of DFA to study the soil geochemistry. The soil geochemical data were treated as compositional data (CoDa). Discriminant function analysis was used in these studies as a robust statistical approach to explore the various factors responsible for different soil geochemical compositions between diverse lithologies, soil types and CLC classes (AITCHISON, 1986). Furthermore, con sistent geochemical maps were created using GIS to visua­ lise the spatial relationships between the models. The main objec tive of the study, based on the DFA results, was to investi­ gate relationships between the groups within these three models, i.e. underlying geology, soil and land use and land cover. 2. MATERIALS AND METHODS 2.1. Study area The Republic of Croatia is geographically located in Central and Southeast Europe, between the latitudes 42° N and 47° N, and longitudes 13° E and 20° E. It covers an area between the central Danube Basin and the central Mediterranean and serves as a bridge between Central and Eastern Europe and the Mediterranean (Fig. 1). Geomorphologically, Croatia can be divided into three different regions, each of which has characteristic climatic and hydrological features: the Pannonian Basin (the northern part), the mountainous region (mostly the central part) and the Adriatic coast (mostly the southern part). The Pannonian Basin features a continental climate with varied atmospheric conditions. FILIPČIĆ (2023) applied the climate regionalisation by KÖPPEN (1918) to Croatia, covering the interval 1961–1990. The Pannonian Basin is assigned to the Cfb region (C – temperate, f – no dry season, b – warm summer). The average temperature in July is 20–22 °C, while the average temperature in January is 0 °C to –3 °C. The amount of precipitation in this area shows a gradient that decreases from the western part of the Pannonian Basin with 1,000–1,100 mm into the central part with 700 to 1,000 mm and to the easternmost part, which is drier with 300 to 700 mm (ZANINOVIĆ et al., 2008). 2.1.1. Geological settings The area of northern and northwestern Croatia geotectonically belongs to the area of the southwestern Pannonian Basin. During the long geological history of northern Croatia (documented over 600 million years), the oldest rocks passed through many orogenic cycles and large proportions of them have been lost (either subducted and/or eroded), so that today we have relatively few surface outcrops of these oldest rocks from the Proterozoic and Palaeozoic (Fig. 2a, b; Table 1). Furthermore, a large part of northern Croatia (about 60%) is covered by the youngest, Quaternary formations (loess, Figure 1. Geographic location of the study area of the Pannonian Basin, Croatia. G eologia C roatica 153Pjanić et al.: Geochemical Map of Croatia – Key findings and legacy of geochemical mapping in the Pannonian Basin Figure 2. a) Geological map (HGI-CGS, 2009); b) the map showing the geological criterion (GEO Code); c) pedological map (BOGUNOVIĆ, 1998); d) the map showing the soil criterion (FAO Code); e) land use and land cover map (CORINE LAND COVER, 2018) and f) the map showing land use and land cover criterion (CLC Code) of the study area in Pannonian Basin, Croatia. The dots in the maps show that GEO, FAO and CLC are also locations of the sampling sites. G eo lo gi a C ro at ic a 154 Geologia Croatica 78/2 aeolian sands, alluvial, lake and marsh sediments) so that older rocks emerge from beneath this young cover as the so­called “Inselgebirge” (VELIĆ & VLAHOVIĆ, 2009). The oldest rocks in that area are Proterozoic (Neopro tero­ zoic) metamorphic rocks, and they are observed on the Psunj and Krndija Mts. The main mass of these rocks consists of various types of paragneisses and orthogneisses, and their source rocks were volcanogenic­sedimentary formations. These rocks were metamorphosed, ranging from chlorite to epidote­amphibolite facies, between 650–550 million years Table 1. Legend of the Geological Map of the Pannonian Basin (from the Geological Map of the Republic of Croatia, scale 1:300,000) (HGI-CGS, 2009). Lithostratigraphic unit Group according to Grouping criteriaNumber of unit Name 58a Deluvial-proluvial deposits (Holocene) GEO12 58b Alluvial deposits (Holocene) GEO13 57b Marsh deposits (Holocene) GEO10 56 Aeolian sands (Holocene) GEO11 54a Terrestrial loess (Pleistocene) GEO10 54b Marsh loess (Pleistocene) GEO10 53a Fluvial deposits (Pleistocene) GEO13 52 Plioquaternary clastic deposits GEO9 50 Paludina beds (Zanclean, Piacenzian) GEO9 49 Sands and clays (Pliocene) GEO9 48 Clastic sediments and coal (Pannonian) GEO9 47 Limestone-clastic sediments (Sarmatian, Pannonian) GEO8 46 Lithothamnium limestone and clastic deposits with volcanics (Badenian) GEO8 45a Magmatic rocks (Karpatian, Badenian): andesites * 45b Magmatic rocks (Karpatian, Badenian): basalts * 44 Clastic and carbonates with clastic (Ottnangian, Karpatian) GEO7 43 Clastic with volcanic sediments (Egerian, Eggenburgian) GEO7 38 Carbonate flysch and clastic (Palaeocene, Eocene) GEO6 37a Magmatic rocks (Upper Cretaceous – Palaeogene): basalts * 37b Magmatic rocks (Upper Cretaceous): rhyolites * 37c Magmatic rocks (Upper Cretaceous): granites * 36 Carbonate clastic (mainly flysch) and “Scagalia” limestones (Upper Cretaceous) GEO6 35 Hemipelagic and turbidite deposits (Lower Cretaceous) GEO6 34 Rudists limestons (Cenomanian, Maastrichtian) GEO6 32 Parametamorphic rocks (Middle Jurassic) * 31a Ophiolite rocks (Middle, Upper Jurassic): ultramafites GEO2 31b Ophiolite rocks (Middle, Upper Jurassic): magmatites GEO2 31c Ophiolite rocks (Middle, Upper Jurassic): sediment rocks GEO2 30 Parametamorphic rocks (Middle Jurassic) GEO2 29 Orthometamorphic rocks (Middle Jurassic) GEO2 28 Calpionellid limestones with cherts (Tithonian, Berriasian) * 27 Platy limestones (Jurassic in general) * 25 Forereef and reef limestones and dolomites (Kimmeridgian, Tithonian) * 20 Dolomites (Upper Norian, Rhaetian) GEO5 17a Magmatic rocks: andesites and basalts (Middle and Upper Triassic) GEO2 16 Clastic and pyroclastic deposits (Middle Triassic) GEO5 15 Carbonate sediments (Middle Triassic) GEO5 14 Induanian and Olenekian deposits (Lower Triassic) GEO4 13b Evaporites and clastic deposits (Upper Permian): clastic GEO4 12 Magmatites (?Permian) * 11 Granites (Permian) * 10 Predominantly clastic deposits (Carboniferous, Permian) GEO3 8 Hercynian semimetamorphic complex (Devonian, Carboniferous, Permian) GEO3 7 Clastic and carbonate sediments (Devonian, Carboniferous) GEO3 6 Parametamorphites (Palaeozoic, Triassic) GEO3 5 Orthometamorphites (Palaeozoic, Triassic) GEO2 4 Granitic rocks (Ordovician, Silurian, Devonian) GEO1 3 Complex of metamorphic rocks (Ordovician, Silurian, Devonian) GEO1 2 Progressive metamorphic series (Ordovician, Silurian, Devonian) GEO1 1 Complex of metamorphic rocks (Precambrian) GEO1 *no sampling sites in these units G eologia C roatica 155Pjanić et al.: Geochemical Map of Croatia – Key findings and legacy of geochemical mapping in the Pannonian Basin ago (JAMIČIĆ & CRNKO, 2009). Since then, there have been no documented rock records in the research area for about another 100 million years. The oldest rocks of the Early Palaeozoic belong to the Or­ dovician. These rocks comprise migmatite gneisses, biotite­ and biotite­muscovite gneisses, mica schists, granitoid rocks, chlorite schists, amphibolites and amphibole schists, pegma­ tites and aplites, peridotites and serpentinites, which are found in the eastern and central parts of Papuk Mt., in the area of Ravna Gora Mt. and the western part of Krndija Mt. The larg­ est part of this complex represents a progressive metamorphic series (PAMIĆ & LANPHERE, 1991), meaning that all tran­ sitions from low to medium metamorphic varieties (chlorite to amphibolite facies) can be found in the terrain. Rocks belonging to the younger Palaeozoic can be found on the Medvednica, Trgovska Gora, Papuk, and Psunj Mts. and to a lesser extent on the Žumberak and Ivanščica Mts. They consist of mostly non­metamorphic to low­metamorphic rocks (slates, slate phyllites, marbles, quartz sericite and quartz chlorite schists). Mesozoic rocks (Triassic and Jurassic) can be observed on almost all the mountains of Northern Croatia (Ivanščica, Kalnik, Medvednica, Žumberak, Trgovska Gora, Psunj and Papuk Mts.). They are composed of sedimentary rocks depos­ ited, initially in continental environments (mainly clastites – Lower Triassic), and after deepening of the sedimentary basin, on a carbonate platform (limestones and dolomites – Middle Triassic to Middle Jurassic). Contemporaneously, in the deeper parts of the ocean space, during the Middle and Upper Trias­ sic and the Lower and Middle Jurassic, rocks of the ophiolitic complex were formed, namely: ultrabasic rocks (peridotites), gabbros, amphibolites, diabases, and pillow lavas, which are interstratified with deep­sea siliceous radiolarian muds. The belt containing these rocks extends from the town of Dvor, through the southwestern parts of Zrinska Gora Mt. and the north­western parts of Trgovska Gora Mt., continuing at a depth of about 800 m below the Vukomeričke Gorice Hills (HGI­CGS, 1992) to the SE part of the Samoborska Gora Mt., and across the NW part of Medvednica Mt. and the SE part of the Ivanščica and Kalnik Mts., continuing towards the north­ east into Hungary (Fig. 2a, b; Table 1). After the closure of the oceanic domain (due to subduc­ tion) from the Late Jurassic to the Early Cretaceous, foreland and piggyback basins were formed, and material eroded from the rising parts of the newly emerging land was deposited in them. Outcrops of these rocks can be found on the Trgovska and Zrinska Gora Mts. in Banovina province, then on the Žumberak, Medvednica, Kalnik and Ivanščica Mts., mostly in the area of the distribution belt of ophiolite rocks. This is fol­ lowed by the intrusion of the Moslavačka Gora granite into the Ordovician granitoid crust, built mainly of muscovite and bio­ tite in different proportions, feldspar and quartz, and belongs to two­mica S­type granite (STARIJAŠ et al., 2010). Transgressive conglomerates, siliciclastic and carbonate sandstones and Scaglia limestones were deposited in the sed­ imentary basin formed in the Upper Cretaceous. Sedimenta­ tion continued with carbonate turbidites until the Palaeocene. In Northern Croatia, during formation of the Paratethys (Eocene–Oligocene), two sedimentary basins developed, the Hrvatsko Zagorje Basin and the Northern Croatian Basin (PAVELIĆ & KOVAČIĆ, 2018), in which sediments were first deposited in terrestrial and freshwater environments, and then, in the Middle Miocene, the sedimentation conditions changed and Lithothamnium limestones were deposited, so that the younger deposits up to the Quaternary gradually change from carbonates to poorly consolidated clastic deposits (siltites, sands, fine­grained conglomerates) (Fig. 2a, b; Table 1). A large part of the central area of northern Croatia is cov­ ered by loess of Pleistocene age (sandy clayey silt, fine­grained sand, silty clay, sandy clay) and marsh loess (clayey sandy silt). The mineral composition indicates its origin from the Alps (MUTIĆ, 1975). The youngest sediments (Holocene) consist of deluvial­proluvial deposits (silts, sands and gravels), marsh sediments (clays, silts, peat) and alluvial deposits (clays, silts, sands and gravels). 2.1.2. Soil The pedological composition of the Pannonian Basin in Croatia is heterogeneous due to considerable differences in the mineral composition of the subsoil, relief and water saturation. The FAO UNESCO pedological map of Croatia contains 32 pedological cartographic units, as well as urban (cities) and aquatic regions (Fig. 2c, d; Table 2; BOGUNOVIĆ et al., 1998). BAŠIĆ et al. (2007) have delineated three primary agricultural regions within the Croatian agrosphere, each distinguished by unique climatic conditions and other characteristics: the Pannonian, the Mountainous and the Adriatic regions, which in turn are subdivided into several subregions. The Pannonian region comprises the eastern, central, western and northwest subregions. The eastern part of the Pannonian region includes Baranja and Eastern Slavonia, located in the easternmost part of Croatia (Fig. 2c, d; Table 2). This area, where loess deposits predominate, favours the formation of very fertile soils such as Chernozem, Cambisol Eutric and Regosol. The Central Pannonian subregion, which includes Western Slavonia, Podravina, Bilogora and Central Posavina, is characterised by the predominance of hydromorphic soils. Stagnosols and Gleysols are the most widespread, followed by Luvisol on loess, Cambisol Eutric, Leptosol on marl and Fluvisols. The western subregion, which lies on the westernmost edge of the Pannonian Basin plain, includes the extensive urban area of Zagreb. The predominant soil types include Luvisol on loess and Stagno­gley, Fluvisol along the Drava River and Gleysols. There are also Eutric Cambisol on loess and Leptosol on marl and soft limestone. Finally, the northwestern subregion comprises the north­west of Croatia, in particular the areas of Zagorje, Varaždin and Međimurje. Due to the heterogeneity of the parent material and the different topography, this subregion has a variety of soil types. In the numerous valleys formed by rivers and streams, Gleysols predominate, followed by Leptosols on marl, Stagnosols and Humo Fluvisols, as well as Regosols and technogenic soils. Technogenic soil is the pedological term used for anthropogenically influenced soils in urban areas (BOGUNOVIĆ et al., 1998). G eo lo gi a C ro at ic a 156 Geologia Croatica 78/2 Table 2. Legend of the FAO UNESCO Pedological map of the Republic of Croatia, scale 1,000,000 of Pannonian Basin (BOGUNOVIĆ et al., 1998). Number of unit Dominant unit Other units Group according to Grouping criteria 1 Mollic and Calcaric Fluvisols (FLm+FLc) Mollic Glaysols (GLm) FL 2 Eutric, Mollic and Calcic Glaysols (GLe+GLm+GLk) Gleyic and Stagnic Podzoluvisols (PDg+PDj) Mollic Fluvisols (FLm) Calcaric Fluvisols (FLc) Calcaric Regosols (RGc) GL 3 Eutric and Calcic Gleysols (GLe+GLk) Fibric Histosols (HSf) GL 4 Mollic Gleysols (GLm) Eutric Gleysols (GLe) Gleyic and Stagnic Podzoluvisols (PDg+PDj) GL 5 Calcaric Regosols (RGc) Lithic Leptosols (LPq) Rendzic Leptosols (LPk) Chromic Cambisols (CMx) Mollic Leptosols (LPm) GL 6 Calcaric Regosols (RGc) Calcic Gleysols (GLk) Calcic Chernozems (CHk) GL 7 Lithic Leptosols (LPq) Mollic Leptosols (LPm) Rendzic Leptosols (LPk) Chromic Cambisols (CMx) LP 8 Rendzic Leptosols (LPk) Aric Anthrosols (ATa) Calcaric Regosols (RGc) Calcic Gleysols (GLk) Stagnic Podzoluvisols (PDj) LP 9 Rendzic and Mollic Leptosols (LPk+LPm) Chromic Cambisols (CMx) Chromic Luvisols (LVx) Aric Anthrosols (ATa) LP 10 Mollic Leptosols (LPm) Chromic Cambisols (CMx) Rendzic Leptosols (LPk) Chromic Luvisols (LVx) LP 11 Umbric and Dystric Leptosols (LPu+LPd) Dystric Cambisols (CMd) Cambic Podzols (PZb) LP 12 Calcaric and Cambic Arenosols (ARc+ARb) Calcic Luvisols (LVk) Stagnic Podzoluvisols (PDj) Dystrict Cambisols (CMd) Mollic and Eutric Leptosols (LPm+LPe) * 13 Eutric and Calcic Vertisols (VRe+VRk) Rendzic Leptosols (LPk) Calcaric Regosols (RGc) Chromic Cambisols (CMx) * 14 Eutric Cambisols (CMe) Calcic Chernozems (CHk) Aric Anthrosols (ATa) CM 15 Eutric Cambisols (CMe) Albic Luvisols (LVa) Mollic Fluvisols (FLm) Eutric Gleysols (GLe) Aric Anthrosols (ATa) CM 16 Eutric and Dystric Cambisols (CMe+CMd) Eutric and Umbric Leptosols (LPe+LPu) Albic Luvisols (LVa) Stagnic Podzoluvisols (PDj) Eutric Gleysols (GLe) Calcaric and Eutric Regosols (RGc+RGe) CM 17 Dystric and Humic Cambisols (CMd+CMu) Cambic and Haplic Podzols (PZb+PZh) Umbric Leptosols (LPu) CM 18 Humic and Dystric Cambisols (CMu+CMd) Umbric Leptosols (LPu) Haplic Lixisols (LXh) CM 19 Chromic Cambisols (CMx) Chromic Luvisols (LVx) Rendzic and Mollic Leptosols (LPk+LPm) Aric Anthrosols (ATa) Eutric Cambisols (CMe) CM 20 Calcic Chernozems (CHk) Eutric Cambisols (CMe) CH 21 Gleyic and Calcic Chernozems (CHg+CHk) Mollic Gleysolls (GLm) CH 22 Albic Luvisols (LVa) Rendzic and Mollic Leptosols (LPk+LPm) LV 23 Albic and Gleyic Luvisols (LVa+LVg) Stagnic Podzoluvisols (PDj) Eutric, Dystrict and Calcaric Cambisols (CMe+CMd+CMc) Stagnic Lixisols (LXj) Cambic and Calcaric Arenosols (ARb+ARc) LV G eologia C roatica 157Pjanić et al.: Geochemical Map of Croatia – Key findings and legacy of geochemical mapping in the Pannonian Basin 2.1.3. CORINE land cover CORINE (COoRdination of INformation on the Environment) is part of the Copernicus Land Monitoring Service (CORINE LAND COVER, 2018). The main technical parameters for Corine Land Cover mapping are the use of a minimum mapping unit of 25 hectares, a minimum width of the linear elements of 100 metres, geometric accuracy of the CLC data of more than 100 m and a thematic accuracy of ≥ 85%. The standard CLC nomenclature comprises 44 land cover classes. These are grouped in a three­level hierarchy. Level 1 consists of five main categories: 1) artificial surfaces, 2) agricultural areas, 3) forests and semi­natural areas, 4) wetlands, and 5) water bodies. At levels 2 and 3, the main categories were subdivided into more detailed classes (Fig. 2e, f; Table 3). 2.2. Sampling and sample preparation Sampling of the topsoil was carried out based on the systematic sampling plan according to ISO 10381­1 (ISO, 2002a) and ISO 10381­2 (ISO, 2002b). The locations of the sampling points and the sampling density throughout the Republic of Croatia were determined by laying out regular square cells with an area of 25 km2 (Fig. 2b, d, f). Five subsamples (0–25 cm depth) were taken from the centres and corners of these cells of ≈400 m² (up to 15% deviation allowed). These were composite samples totalling 3–5 kg stored in PVC bags. The samples were air dried, sieved to a fraction of < 0.063 mm and homogenised in an agate mortar. According to the protocol (HALAMIĆ et al., 2000), the sampling documentation included field observations and administrative details (sample number, soil type, coordinates, altitude, images of the surroundings, date of sampling) as well as field observations (vegetation, contaminants, relief, depth of soil profile, colour, texture, subsoil lithology), which contribute to the database of the digital geochemical atlas of Croatia. 2.3. Laboratory analyses and quality control The homogenized soil samples were dissolved in a mixture of four acids HF­HCl­HNO3­HClO4. The solutions were analysed by Inductively Coupled Plasma Mass Spectroscopy using a Perkin Elmer Elan 6000 or 9000 ACME Analytical Laboratories (now Bureau Veritas Mineral Laboratories), Vancouver, Canada (ACME, 2007) for a range of 41 elements (Ag, Al, As, Au, Ba, Be, Bi, Ca, Cd, Ce, Co, Cr, Cu, Fe, Hf, K, La, Li, Mg, Mn, Mo, Na, Nb, Ni, P, Pb, Rb, S, Sb, Sc, Sn, Sr, Ta, Th, Ti, U, V, W, Y, Zn, and Zr). Mercury analysis was performed using aqua regia extraction by flameless atomic adsorption spectrometry (FAAS). Analytical precision was checked by duplicate analyses of both certified reference samples and randomly selected soil samples (every 20th sample in the batch), resulting in an average coefficient of variation of approximately 5%. Accuracy was checked with certified geological reference materials, i.e., GXR­2 and GXR­5 (ACME Labs.), and SJS­1 (USGS). The accuracy of most elements analysed in the reference soils was within ± 10% of the certified values (HALAMIĆ & MIKO, 2009). 2.4. Data processing A geochemical database of the geochemical, geological, pedological and other relevant data was created in the GIS software ESRI® ArcGIS™ 10.2.1. Number of unit Dominant unit Other units Group according to Grouping criteria 24 Chromic Luvisols (LVx) Chromic Cambisols (CMx) Mollic Leptosols (LPm) Dystrict Cambisols (CMd) Albic Lixisols (LXa) LV 25 Stagnic and Gleyic Podzoluvisols (PDj+PDg) Eutric Cambisols (CMe) Albic and Stagnic Luvisols (LVa+LVj) Eutric and Mollic Gleysols (GLe+FLm) Dystrict Cambisols (CMd) PZ 26 Cambic Podzols (PZb) Humic Cambisols (CMu) Haplic Podzols (PZh) Umbric Leptosols (LPu) PZ 27 Haplic and Cambic Podzols (PZh+PZb) Humic Cambisols (CMu) Umbric Leptosols (LPu) PZ 28 Humic Acrisols (ACu) Chromic Cambisols (CMx) Rendzic Leptosols (LPk) Dystrict Cambisols (CMd) * 29 Fibric and Terric Histosols (HSf+HSs) Calcic and Mollic Gleysols (GLk+GLm) * 30 Aric Anthrosols (ATa) Cumulic Anthrosols (ATc) URB 31 Mollic Gleysols (GLm) Gleyic Solonchaks (SCg) GL 32 Gleyic and Albic Luvisols (LVg+LVa) Stagnic and Glayic Solonetz (SNj+SNg) Calcic Solonetz (SNk) * 33 Urban areas – URB 34 Water surfaces (lakes, ponds, sea) – *no sampling sites in this unit G eo lo gi a C ro at ic a 158 Geologia Croatica 78/2 2.4.1. Statistical processing 2.4.1.1. Compositional data and log-ratio analysis The analysed dataset is composed of 28 variables exploited consistently in various earlier investigations pursuing the low­ density soil sampling during the geochemical baseline mapping in Croatia (e.g., MIKO et al., 2001; PEH et al., 2010; HALAMIĆ et al., 2012). It includes Al, As, Ba, Ca, Co, Cr, Cu, Fe, Hg, K, La, Mg, Mn, Na, Nb, Ni, P, Pb, Sc, Sn, Sr, Th, Ti, V, Y, Zn, Zr, together with pH as the measure of soil acidity. This element suite represents the input data (explanatory variables) in discriminant function analysis used to study the impact of various soil­forming factors in the PAN part of Croatia. Descriptive statistics for the entire dataset (Min, Max, Mean and Median) are shown in Table 4 indicating, albeit, information that is proper for comparison purposes only, because the data displayed therein represent relative instead of absolute values. Namely, a well­known fact that soil geochemical data represent a summary of compositional data categorically interferes with their use in the raw (compositional) form in any statistical analysis (EGOZCUE & PAWLOWSKY­ GLAHN, 2006). By its very nature, CoDa involves the mathematical particularity that all variables (parts) in each case (analysed sample) must be positive, confined to a constant sum a priori specified as 100%, 106 mg/kg, or 1.0. On account of this unit­sum constraint, CoDa can be genuinely displayed only in a limited sample space (compositional space) also named simplex, composed of D parts or components (geochemical variables). Thus, a suite of D­part composition (SD) fills a restricted part (between null and 100%) of a D­dimensional real space (RD), building a subset of its vectors (PAWLOWSKY­GLAHN & EGOZCUE, 2006; BUCCIANTI, 2013; BUCCIANTI & GRUNSKY, 2014). The philosophy of the simplex as the regular sample space to inhabit CoDa is articulated with the following algorithm (PAWLOWSKY­ GLAHN et al., 2007; BUCCIANTI, 2013): SD = x x x x x x i D xD i i D i= ( ) > ( = ) ={ }=∑1 2 3 1 0 1 2 3, , ,..., : , , ,..., , (1) where x is a constraint­sum constant; x1, x2, x3, . . ., xD are components of the composition x; and 1, 2, 3, . . . , D are parts of the composition x. Simplex can be “unfolded” into the structure of the Euclidean vector space only after a suitable transformation of its components. One of the several transformations normally exploited in CoDa analysis is the centred log ratio (clr) of the raw (compositional) data. Originally introduced by AITCHISON (1986), this transformation is extensively utilised in multivariate statistics, including this study. The application of clr coefficients is considered of crucial importance in multivariate analyses such as DFA, since it maintains the Table 3. Legend of Corine Land Cover 2012 (CORINE LAND COVER, 2018) raster data of the Pannonian Basin (European Environment Agency). Level 1 CLC2018 Description Level 2 CLC2018 Description Group according to Grouping criteria 1 Artificial surfaces 111 Continuous urban fabric * 112 Discontinuous urban fabric ARTS 121 Industrial or commercial units ARTS 122 Road and rail networks and associated land * 123 Port areas * 124 Airports * 131 Mineral extraction sites * 132 Dump sites * 133 Construction sites * 141 Green urban areas * 142 Sport and leisure facilities ARTS 2 Agricultural areas 211 Non-irrigated arable land AGRS 212 Permanently irrigated land AGRS 221 Vineyards AGRS 222 Fruit trees and berry plantations AGRS 231 Pastures AGRS 242 Complex cultivation patterns AGRS 243 Land principally occupied by agriculture, with significant areas of natural vegetation AGRS 3 Forest and seminatural areas 311 Broad-leaved forest FSNA 312 Coniferous forest FSNA 313 Mixed forest FSNA 321 Natural grasslands FSNA 324 Transitional woodland-shrub FSNA 331 Beaches, dunes, sands * 4 Wetlands 411 Inland marshes WETL 5 Water bodies 511 Watercourses WETL 512 Water bodies WETL *no sampling sites in this unit G eologia C roatica 159Pjanić et al.: Geochemical Map of Croatia – Key findings and legacy of geochemical mapping in the Pannonian Basin original distances between the corresponding components (EGOZCUE & PAWLOWSKY­GLAHN, 2006; TOLOSANA­ DELGADO, 2012). The common problem of singularity, intrinsic to the clr­transformed covariance matrix, can be easily circumvented because DFA operates with the reduced data matrix, not contingent upon the full rank covariance matrix (DAUNIS­I­ESTADELLA et al., 2011). Actually, by the scale invariance property of the log­ratio approach, more coherent results are obtained in the case when only a subset of the variables (sub­composition) instead of full composition is used in the analysis (MERT et al., 2018), a procedure involving exclusion of one composition (variable) after transformation. Given that clr­transformed data assume the role of unconstrained actual vectors in the real (Euclidean) space, the Mahalanobis distances (MD) between cases and groups remain invariant irrespective of which component is removed from the analysis (BARCELÓ­VIDAL et al., 1999). The clr­coefficients can be calculated from the following expression: clr x x g x x g x x g x x g x D( ) ln ( ) , ln ( ) , ln ( ) ,..., ln ( ) =     1 2 3 (2) where x1, x2, x3, . . ., xD are parts (compositions) of x, while g(x) represents their geometric mean. Note that the clr coefficients are computed by dividing each component (such as the geochemical percentage, or mg/ kg) by the geometric mean of all components involved in the analysis and obtaining their logarithm in the last analysis (AITCHISON, 1986), as displayed in equation (2). The clr­ transformed variables represent dimensionless numbers (ratios) unsuitable for direct cross­comparison, functioning simply as input data for DFA. The quality explained above adheres to the logic of the log­ratio approach itself, substantiating further communication about the algebraic­geometric structure of the sample space. This is important for a better understanding of the results of the DFA presented in this study, especially in explaining why separate (variable and sample) scatterplots are used instead of biplots. As elucidated by AITCHISON & EGOZCUE (2005) and described in GALOVIĆ & PEH (2016), compositional parts have a twofold nature: they can be portrayed both as raw compositional data such as percentages (or mg/kg) of geochemical compositions, using vectors of parts and as coordinates (scalars) in the Cartesian (orthonormal) coordinate system with Euclidean metrics. In the log­ratio approach, they are not considered as simple conversions of the original data into their logarithms for data normalization, but as coordinates. Seen in this light, straightforward transposition of the original data (raw compositions) to coordinate scatterplots or biplots can easily produce erroneous interpretations. One of the most tenacious fallacies that sometimes favours the application of classical (non­compositional) statistical methods Table 4. Descriptive statistics for raw (compositional) geochemical data. Min Q1 Med Q3 Max MAD g Al (%) 1.46 5.79 6.41 7.20 11.81 0.71 6.45 As (mg/kg) 0.50 6.00 9.00 12.00 92.00 3.00 8.10 Ba (mg/kg) 58.00 366.00 403.00 451.00 3300.00 40.00 406.33 Ca (%) 0.07 0.48 0.72 1.41 26.79 0.31 0.94 Co (mg/kg) 3.00 8.00 10.00 13.00 36.00 2.00 10.56 Cr (mg/kg) 28.00 63.00 75.00 92.00 524.00 13.35 77.48 Cu (mg/kg) 3.00 14.10 19.05 28.20 248.00 6.05 19.89 Fe (%) 0.60 2.60 2.99 3.56 6.81 0.46 3.05 Hg (μg/kg) 5.00 30.00 40.00 65.00 4535.00 15.00 44.51 K (%) 0.33 1.49 1.65 1.82 3.42 0.17 1.64 La (mg/kg) 9.00 33.00 38.00 42.30 71.20 4.60 37.12 Mg (%) 0.23 0.61 0.75 1.07 7.52 0.18 0.86 Mn (mg/kg) 131.00 462.00 579.00 730.00 5619.00 131.00 579.60 Na (%) 0.11 0.76 1.03 1.24 3.21 0.24 0.94 Nb (mg/kg) 1.00 8.00 9.00 10.40 24.10 1.00 9.13 Ni (mg/kg) 9.20 25.90 32.60 45.00 427.00 8.60 34.81 P (%) 0.02 0.05 0.07 0.09 0.41 0.02 0.07 Pb (mg/kg) 14.00 22.00 25.75 32.00 699.00 4.75 28.11 Sc (mg/kg) 2.00 9.00 10.00 12.00 25.00 2.00 10.09 Sn (mg/kg) 0.10 3.30 4.00 5.00 70.00 1.00 3.99 Sr (mg/kg) 35.00 100.00 111.00 127.00 1090.00 12.00 115.38 Th (mg/kg) 2.00 10.00 11.60 12.90 20.00 1.40 11.15 Ti (%) 0.08 0.35 0.40 0.43 1.13 0.04 0.39 V (mg/kg) 22.00 79.00 89.00 105.00 238.00 12.00 91.58 Y (mg/kg) 3.00 13.00 16.60 19.60 48.00 3.40 15.56 Zn (mg/kg) 28.00 62.00 72.00 93.00 1432.00 14.00 78.92 Zr (mg/kg) 9.00 26.00 36.00 42.30 1583.00 7.66 33.20 pH 2.03 5.01 5.72 6.81 8.09 0.86 5.79 Note: Q1, Med, and Q3 are the sample quartiles (25th, 50th and 75th percentile); MAD is median absolute deviation; g is geometric mean G eo lo gi a C ro at ic a 160 Geologia Croatica 78/2 on CoDa arises from the fact that the latter sometimes simply reaffirms what has already been confirmed by traditional methods. In this case, though, attention is required because “... either we have been lucky with our traditional methods, or at least the new methodology must be correct in this case” (AITCHISON, 2008). 2.4.1.2. Discriminant function analysis – creating a strategy DFA is one of the traditional multivariate statistical techniques, which is extremely helpful in creating the predictive model of a two­ or multiple­group discrimination based on a suite of independent (predictor) variables. In its various practical considerations in geosciences, DFA is commonly used when geological rationale requires the application of some autonomous discriminating criterion concerning the variables in the analysed suite of data. It is an efficient tool in handling the large number of quantitative attributes, reducing problems with organization, distinction, or comparison of sizeable databases, to gain a higher degree of understanding of the various underlying geological, geochemical or environmental controls. In addition, the data managed in this way may also result in mapping properties that help to clarify the inherent relationships among the original variables. The aims and principles of DFA are described in many statistical textbooks (e.g. DAVIS, 1973, 1986; DILLON & GOLDSTEIN, 1984; ROCK, 1988; REIMANN et al., 2008) and are regularly put forward by the present authors in various geochemical, environmental, sedimentological, bauxite, and other studies (PEH et al., 2008; PEH & HALAMIĆ, 2010; KOVAČEVIĆ GALOVIĆ et al., 2012; HALAMIĆ et al., 2012; PEH & KOVAČEVIĆ GALOVIĆ, 2014, 2016; GALOVIĆ & PEH, 2016; GRIZELJ et al., 2017; ŠORŠA et al., 2018a; BRUNOVIĆ et al., 2019; HASAN et al., 2020; GIZDAVEC et al., 2022; ILIJANIĆ et al., 2023; GALOVIĆ et al., 2024). As a multivariate method, it was applied in this study to process a vast body of data (1,254 samples) in a manner accommodating the most efficient concordance between the soil geochemical signature and various aspects of the surrounding soil environment in the Pannonian part of Croatia. Characterisation of the grouping criteria is critical in this regard because geochemical patterns in the sampling media (upper soils in this case), typically follow the overall perspective on a regional scale – geological, environmental and other intrinsic controls predominating in the investigated area (Croatian Pannonian region). These principles are autonomous concerning the analysed (independent) variables (see, e.g., ROCK, 1988). A discriminant analysis panel from the STATISTICA statistical software package, version 7.1 (STATSOFT INC., 2006), was used to achieve the best possible separation between the á priori defined groups and to determine the geological/ecological/ pedological factors responsible for the structure of the input data. In this procedure, the applied statistical package (linear discriminant function) automatically finds the best variable patterns, allowing computed discriminant functions to make the maximum contribution to the multi­group discrimination. In this way, casual confusion and possible misinterpretation of group memberships in the computed model is prevented, at least to the extent that DFA is applied for confirmatory purposes. One of the most decisive standards appropriate to the group description in this case is represented by the underlying geology (lithology). It draws heavily on earlier research work (HASAN et al., 2020) that built the strong geochemical contrast between the soil geochemistry of the two geological/ structural regions in Croatia roughly delineated as DIN and PAN. In particular, that study demonstrated that even the lithologically monotonous terrain, such as the Dinaric­coastal region, which is characterized by predominantly carbonate bedrock, appears sufficiently diverse to affect the geochemical signal in the overlying soils. Furthermore, it showed that recent investigations concerning the GEMAS Project (Geochemical mapping of agricultural and grazing land soil) (REIMANN et al., 2014) validated the suitability of the land cover classes, borrowed from the CLC inventory, in the search for environmental impacts on the geochemical composition of soils. Finally, soils by their very nature are a repository of various processes, past or present, participating in soil formation and leaving distinctive geochemical signals at both local and regional levels in the final analysis. To this extent, various soil classifications, such as FAO (developed by the Food and Agriculture Organization of the United Nations), offer valuable generalizations about pedogenesis with respect to the interaction between the key soil­forming factors. Following the guidelines declared above, three substantive issues in this study were delineated with respect to the grouping strategy: GEOLOGY, CLC, and SOIL. The first two of these include a different number of classes depending on the nature of the grouping variables which stem, at least partly, from the familiar ‘clorpt’ state equation: S = f (DFcl, o, r, p, t); published by JENNY (1941), while the last one represents the final solution of the entire process. In other words, the soil is the function of climate, organisms, relief, parent material and time, explaining the role of variables (state factors) in the process of soil formation (e.g., PHILLIPS, 1998, 2002; BOCKHEIM et al., 2005). Later on, this state­factor model (unsolvable in its original form) was extended to include the ecosystem, soil, vegetation and fauna (e.g., BOCKHEIM & GENNADIYEV, 2010) and finally reviewed in the work on soil complexity and pedogenesis (PHILLIPS, 2017). In this study, the groups are formed using the following references: 1) the GEOLOGY model is based on the general geology of the investigated area introduced from the Geological Map of the Republic of Croatia, 1:300,000 (HGI­CGS, 2009) (13 groups); 2) the CLC model exploited the most general level of standard CLC classification (Label 1) from the CLC Project 2018 (4 groups; where the CLC2018 groups labelled as 4 and 5 were joined into one), while; 3) the SOIL model made use of the FAO UNESCO Pedological map of the Republic of Croatia, scale 1,000,000 of Pannonian Basin (BOGUNOVIĆ et al., 1998; 8 groups). All three models comprise 1,254 valid objects altogether (N) and the same suite of variables (p = 28), as described in Table 5. 2.4.2. The map generation The geochemical maps were created with the ArcGIS extension Geostatistical Analyst. Inverse Distance Weighting (IDW) was used to visualise the maps with the discriminant scores. The method parameters were the same as those used to create the geochemical maps of the elements in the Geochemical Atlas of Croatia (ŠORŠA, 2009). The calculation is performed https://en.wikipedia.org/wiki/Food_and_Agriculture_Organization_of_the_United_Nations https://en.wikipedia.org/wiki/Food_and_Agriculture_Organization_of_the_United_Nations G eologia C roatica 161Pjanić et al.: Geochemical Map of Croatia – Key findings and legacy of geochemical mapping in the Pannonian Basin with a maximum of 9 neighbouring points and a minimum of 6 points. The calculated area is a circle with a diameter of 10 km, choosing a power value of p = 2 to achieve an optimal and balanced mutual influence between the sample points. The results of the discriminant scores are categorised into 8 percentile classes, with the boundaries of each class at the 5th, 10th, 25th, 50th, 75th, 90th and 98th percentile. The spatial distribution of each class is shown on the map in the form of a polygon. To achieve a smoother transition on the maps, a new polygon was used for every 2nd percentile. The colours range from shades of blue for low concentrations to green­yellow­ orange and red for high concentrations. 2.5. Description of the studied groups 2.5.1. Lithological groups Due to the great lithological diversity, and the relatively small distribution of individual types of rocks, we tried to classify them into groups with common characteristics. The idea was to define each group by its lithological similarities and a similar time of formation. In addition, some groups were formed, based on the genetic characteristics of certain types of rocks. For example, group two includes ophiolitic rocks, but also other types of rocks that are genetically related to them (Fig. 2a, b; Tables 1, 4). Group 1 (GEO1) – includes the oldest metamorphic rocks in northern Croatia (Neoproterozoic), which are found in the mountains of the Psunj, Krndija, and a small part of Moslavačka Gora. These are metamorphites from the progressive­ metamorphic series such as ortho­ and paragneisses, amphibolites, amphibole schists, metagabbros, marbles, granite and granitoid deposits, chlorite schists and graphite schists. Group 2 (GEO2) – consists mainly of ophiolitic rocks (Triassic and Jurassic) and sedimentary rocks that are genetically related to them and metamorphic products created by the metamorphosis of ophiolitic rocks. The rocks of this group are distributed in a relatively narrow ophiolite belt that extends from the central Dinarides of Bosnia and Herzegovina over the Zrinska and Samoborska Gora Mts., then over Medvednica and Ivanščica to the Kalnik Mts, where it sinks to the northeast beneath Tertiary and Quaternary sediments and continues into Hungary. Ultramafic rocks (lercolites and serpentinised peridotites), as well as mafic rocks (gabbros, diabases and basalts), are included in this group. Additionally, the group comprises amphibolites and various varieties of green ortho­schists derived from ophiolitic rocks. Associated with this group of rocks are the sedimentary rocks of the ophiolitic mélange (matrix­supported conglomerates, sandstones and siltites) as well as low grade metamorphic rocks including slates, phyllites, metasandstones and quartz schists. Group 3 (GEO3) – contains metamorphic, low grade metamorphic and sedimentary rocks of the Palaeozoic (Devonian and Carboniferous), including Lower Permian deposits (mainly products of the Variscan orogeny). These rocks form parts of the Trgovska Gora, Petrova Gora, Samoborska Gora, Medvednica, Papuk, Psunj and Krndija Mts. The rocks are represented by different types of schists, phyllites, metagreywackes, marbleized limestones, marbles, recrystallized dolomites, shales, siltites, sandstones, fine­ grained conglomerates, argillaceous limestones and dolomites. Group 4 (GEO4) – consists mainly of clastic deposits of the Upper Permian and Lower Triassic. Part of these sediments are products of terrestrial depositional environments. The rocks extend over the Papuk, Krndija, Petrova Gora, Žumberak, Samoborska Gora, Medvednica, Ivanščica, Strahinjščica and Ravna Gora Mts. They comprise shales, siltites, sandstones, conglomerates, marls and biocalcarenites. Group 5 (GEO5) – the rocks formed on the Mesozoic carbonate platform are included in this group. These are predominantly carbonate deposits of Middle Triassic to the Middle Jurassic age. They were detected on the Trgovska and Petrova Gora, Žumberak, Medvednica, Ivanščica, Strahinjščica, Kuna and Ravna Gora, Kalnik, Psunj, Papuk, Požeška Gora and Krndija Mts, and they consist predominantly of limestones and dolomites. In addition, shales, siltites, sandstones, deep­water limestones and tuffs of volcano­ sedimentary origin (Pietra verde) are also represented. Table 5. Grouping criteria. N GEOLOGY n SOIL n CLC n 1 Metamorphic rocks (GEO1) 32 Chernozem (CH) 31 Agricultural surfaces (AGRS) 662 2 Ophiolites (GEO2) 9 Cambisols (CM) 182 Artificial surfaces (ARTS) 49 3 Clastic and carbonate deposits (GEO3) 22 Fluvisols (FL) 145 Forests and seminatural areas (FSNA) 505 4 Induanian and Olenekian deposits (GEO4) 7 Gley (GL) 293 Wetlands (WETL) 38 5 Carbonate sediments (GEO5) 16 Leptosols (LP) 83 6 Carbonate flysch and clastic sediments (GEO6) 12 Luvisols (LV) 257 7 Clastic and carbonate deposits with volcanics (GEO7) 32 Podzols (PZ) 250 8 Limestone and clastic deposits with volcanics (GEO8) 77 Urban soils (URB) 13 9 Clastic deposits, sands, clays, coal (GEO9) 156 10 Loess and marsh deposits (GEO10) 489 11 Aeolian sands (GEO11) 32 12 Deluvial-proluvial deposits (GEO12) 38 13 Alluvial deposits (GEO13) 332 Total 1254 1254 1254 Note: N = number of groups; n = number of cases in each respective group (∑N = 1,254) G eo lo gi a C ro at ic a 162 Geologia Croatica 78/2 Group 6 (GEO6) – Sediments from the Lower Cretaceous to the Palaeogene are included in this group and are found on the Zrinska Gora, Žumberak, Medvednica and Požeška Gora Mts. These deposits consist of turbidites (rhythmic alternation of shales, siltites, sandstones and conglomerates), then greywacke sandstones, calcarenites, calcrudites, conglomerates and cherts. In addition, hemipelagic limestones of the Scaglia facies and shallow marine biogenic limestones also appear. Group 7 (GEO7) – includes clastic sediments and volcanics of Oligocene and Lower Miocene ages, which are distributed on the Kalnik, Ivanščica, Strahinjščica and Kuna Gora Mts. In addition, clastics with carbonate rocks from the upper part of the Lower Miocene were also included in this group. They are found around almost all the mountains in northwestern Croatia. These rocks are made up of sands, sandstones, gravels, conglomerates, marls, clays, siltites, coals and tuffs as well as Congeria limestones. Group 8 (GEO8) – Limestones, clastic deposits with volcanics and carbonate­clastic deposits of the Middle and Lower part of the Upper Miocene belong to this group. They consist of Lithothamnium limestones, argillaceous limestones, micritic and granular limestones, sandstones, conglomerates and marls, rhyolites and acid tuffs. All these rocks surround most of the northern Croatian mountains, and on Dilj Mt. they form its core as well. Furthermore, they were mapped on the Zrinska Gora and Bansko Brdo Mts. in Baranja. Group 9 (GEO9) – comprises sedimentary rocks from the Upper Miocene to the Plio­Quaternary. These are clastic sediments such as marls, siltstones, multi­coloured sands and sandstones, calcite clays, calcarenites, gravels, coal and carbonaceous clays. These deposits are found around all the mountains of northern Croatia, and they fill the Drava and Sava basins and the Karlovac and Bjelovar depressions. Group 10 (GEO10) – this group consists of deposits of terrestrial and marsh loess, which originate from the Alps. They are found on Bilogora Mt., in Moslavina, in the Bjelovar and Ilova depressions, in the Karlovac and Požega valleys, on the Đakovo–Vinkovci–Vukovar plateau, and on the Erdut and Bansko Brdo Hills. Terrestrial loess consists of silts and fine­ grained sands, while marsh loess consists of clayey silts and sands, and sands. Marsh loess was formed by the deposition of aeolian sediment directly into the lake or marsh environment. Group 11 (GEO11) – consists only of aeolian sands, which were formed by the blowing of fine­grained river (terrace) sediments of the Mura and Drava Rivers by the north­ west winds, which explains why these deposits are found only along the southern edges of the Drava and Mura basins. These sediments are spread in the northern part of Međimurje province from Mursko Središće in the northwest to the Goričan towns in the southeast. South of the Drava River, these deposits are found in a relatively narrow belt from the towns of Varaždin to Ludbreg, continuing south­east along the northern slopes of Bilogora Mt. to the Virovitica town, from where they extend eastwards directly along the Drava River to Valpovo. In that area, hectometre­ to kilometre­long dunes with an approximate northwest­southeast orientation are also visible. Aeolian sediment consists of sands and siltstones. Group 12 (GEO12) – includes deluvial­proluvial deposits found in larger areas at the foothills of the Bilogora, Krndija, Moslavačka Gora and Psunj Mts., Erdutsko Brdo Hill and in the area between the towns of Sunja and Dubica in the Banovina province. The composition of these sediments depends directly on the eroded substrate, and they are mostly composed of blocks, rock fragments, semi­rounded gravels, sands and silts. Group 13 (GEO13) – consists of the youngest Holocene sediments, which are still formed today in valleys (floodplains) and riverbeds, and streams in the entire area of northern Croatia. These sediments occupy the largest areas in the valleys of the Sava, Drava, Mura, Krapina, Kupa, Lonja, Česma, Ilova, Pakra, Orljava, Karašica and Danube Rivers. The sediments consist of unconsolidated gravels, sands, silts and clays. 2.5.2. Soil groups BOGUNOVIĆ et al. (1998) harmonised the Croatian pedological nomenclature with the FAO nomenclature and published the FAO UNESCO pedological map of the Republic of Croatia (Fig. 2c; Table 2). Since the pedological data in the European Soil Database (PANAGOS, 2006) are more general than those of BOGUNOVIĆ et al. (1998), the grouping criterion of the latter was used. The criterion for a particular group was the soil type of the Dominant units, except for three soil types that were present in an insignificant number of our sampling sites. These soil types were linked by Other units (1), Regosol was linked to Glaysols, (2) Calcic Chernozem was linked to Cambisol, and (3) Gleyic and Calcic Chernozem were linked to Gleysols (Fig. 2c, d; Tables 2, 4). The unmapped urban areas, such as the city of Zagreb, were identified by the authors as areas with Urban soils. Group FL – Fluvisols from cartographic unit 1 are soils formed by the deposition of soil particles on the floodplains of rivers. The relatively young Quaternary age of the sediments and the repeated flooding and deposition of new particles have prevented pedogenic processes, so Fluvisols are classified as underdeveloped hydromorphic soils. These soils are developed in the alluvial and lacustrine deposits of the Sava, Drava, Mura and Danube Rivers in the Pannonian Basin. Group GL – The Gleysols identified in units 2, 3, 4 and 31 and the Regosols from unit 6 are water­saturated most of the year. Gleysols are mainly found along the rivers such as the Sava, Drava, Mura, Lonja, Kupa, Krapina, Česma, Ilova, Pakra, Orljava and Danube, predominantly in lowland areas. Their water­saturated condition facilitates the formation of gleys, which are characterised by a grey­blue colour caused by the reduction of iron under water­saturated conditions. Regosols, which have been added to this FL group, are shallow, poorly developed soils on loose sediments such as loess, sand, marl and flysch. They are found on loess in Eastern Slavonia (Baranja, Fruška Gora Mt.). Group LP – Leptosols from units 8, 9 and 11 are shallow soils with limited pedogenic development, typically overlying hard rock. These soils show minimal horizon development due to the limited material available for soil formation and are often found in mountainous or rocky terrain. In the Pannonian G eologia C roatica 163Pjanić et al.: Geochemical Map of Croatia – Key findings and legacy of geochemical mapping in the Pannonian Basin Basin, they develop on the slopes of “Inselgebirge” in Slavonija (Dilj, Krdnija, Moslavačka Gora Mts.) and Banovina (Zrinska Gora Mt.), then in the hilly Hrvatsko zagorje and Međimurje provinces. Group CM – Cambisols occur in units 14, 15, 16, 18 and 19, and are moderately developed, young soils. These soils are characteristic of regions where weathering and biological activity are strong but have not yet reached advanced stages of soil maturity. Eutric cambisols are mainly found on the loess plateaus of the Đakovo, Vinkovci, Vukovar and at the “Inselgebirge” (Bansko Brdo and Erdutsko Brdo Hills, and the Dilj, Papuk, Bilogora, Moslavačka Gora, Žumberak Mts). Group CH – Chernozems dominate in cartographic units 20 and 21. They are rich in organic matter and secondary calcium carbonate concretions. These fertile, deep black soils, which occur mainly in eastern Croatia, are known in areas such as Dalj and Zmajevac (GALOVIĆ et al., 2023; POCH et al., 2024). Group LV – Luvisols are found in cartographic units 22 and 23. They are the most common soil type in Croatia, particularly widespread in the Pannonian Basin. A Luvisol is an automorphic soil type with a subsurface horizon with high clay accumulation and high base saturation. Luvisols are particularly common in Međimurje, Hrvatsko zagorje and in the hilly areas of Vukomeričke Gorice and Petrova Gora. These soils also extend across the Ivanščica, Kalnik, Bilogora, Moslavačka Gora, Papuk and Krndija Mts. and are occasionally found in Srijem and central Baranja. Group PZ – Podzols (all sampling sites are located in Podzoluvisols), assigned to cartographic unit 25, are acidic and waterlogged soils in which organic matter accumulates faster than it decomposes, leading to the leaching of minerals and the development of a spodic horizon. These soils occur almost the entire studied area, except in Međimurje and Hrvatsko zagorje provinces, where Podzoluvisols occur only sporadically. Group URB – urban soils; includes soils in urban areas that are not mapped and are labelled as “cities” on the FAO UNESCO pedological maps (BOGUNOVIĆ et al., 1998). Urban soils have developed under anthropogenic influences and could vary considerably in their properties due to disturbances such as construction, landscaping and pollution, which affect their natural soil profile and function. The absence of the cartographic soil units 5, 7, 10, 12, 13, 17, 24, 26, 27, 28, 29, and 32 in the PAN region illustrates the regional specificities in soil formation and distribution, underlining the importance of local geological, hydrological and climatic conditions in pedogenesis. 2.5.3. Land use and land cover groups For the CLC classification, the most general level 1 of the standard CLC classification from the CORINE LAND COVER (2018) raster data of the European Environment Agency was used and summarised in 4 groups. The first three groups correspond to the first three CLC classes: (1) artificial surfaces/urban or built­up areas (ARTS), (2) agricultural land (AGRS) and (3) forests/forest land semi­natural area (FSNA). The fourth group was formed by combining the two CLC Level 1 classes Wetlands and Waterbodies into one group (4) Wetlands (WETL) (Fig. 2e, f; Tables 3, 4). Group ARTS – includes areas with strong anthropogenic influences. The sampling sites were located in the three CLC2018 areas: discontinuous urban fabric (buildings, with associated land, access roads, and car parks), industrial or commercial units (manufacturing, trade, services, and transport, with associated land) and sports and leisure facilities (parks, sports facilities) (Fig. 2e, f; Tables 3, 4). Group AGRS – includes areas with various types of agricultural production: non­irrigated arable land, permanently irrigated land, vineyards, fruit trees and berry plantations, pastures, complex cultivation patterns and land, mainly agricultural land with significant areas of natural vegetation. Group FSNA – represents ground covered with forest and semi­natural areas. This group includes broad­leaved forests, coniferous forests, mixed forests, natural grasslands and transitional woodland shrubs. Group WETL – integrate inland marshes flooded for most of the year, natural and artificial watercourses that serve as drainage channels, and water bodies such as lakes, ponds, reservoirs and canals. 3. RESULTS AND DISCUSSION Results of the DFA are presented in the common table (Table 6) comprising all three exploratory models (tests of residual roots). Before that, the overall significance of their discrimination potential is verified by the appropriate multivariate tests (Table 7) exposing the vanishingly small associated probabilities (p < 0.000 for overall significance of discrimination), a condition required to securely proceed with calculating discriminant functions (DFs). According to various numbers of groups (K) in each model, the total number of DFs is K­1. Notwithstanding the statistical significance of the variation between the observed groups (p­level in Table 6), which defines the dimensionality of the discriminant space, only a smaller number of DFs is handled to explain the natural variation between groups. Natural variation (geological meaning) hidden behind the original data served as a helpful criterion decreasing the number of discriminant axes to only two or three. However, unlike the karst region, where the lithological contrast between carbonate and siliciclastic facies was clear enough to reduce the number of groups to only a few, the Pannonian Basin appears as a true mosaic of lithological diversity. As a result, the first three DFs scarcely exceed half of the total variance comprised in the thirteen observed GEOLOGY groups. 3.1. Functional models – labelling the discriminant functions The technique of labelling DFs is explained elsewhere, including the reasons for using scatterplots instead of biplots in the CoDa analysis (e.g., PEH & KOVAČEVIĆ GALOVIĆ, 2014; GALOVIĆ & PEH, 2016; ŠORŠA et al., 2018a). In essence, the structural (mathematical) model is converted into a functional (process) one, which is geochemical at its core. In G eo lo gi a C ro at ic a 164 Geologia Croatica 78/2 the process, the group centroids (means) are utilised as the proxy for the cloud of single objects in the scatterplot construction to simplify the explanation of the models. Scatterplots of group centroids and variable loadings are built for all three discriminant models using the first two DFs explaining the greatest proportion of the between­group variance. Models are contrasted using the scatter diagrams of the DF1 and DF2 pairs of discriminant functions (orthogonal axes) (Fig. 3). 3.1.1. GEOLOGY model In the GEOLOGY model, generally concerned with the parent­ material state variable (“p”), from the aforementioned “clorpt” state equation, the first discriminant function DF1 explaining one­third of the total variance (33.31%) highlights the youngest (Quaternary) sediments (Table 1; GEO10 – GEO13) in contrast to all the remaining lithologies, a clearly bipolar relationship (Fig. 2a, b; Fig. 3a, b; Fig. 4a, b; Table 6). Undoubtedly, soils developed on alluvium and other Quaternary deposits stand out as a repository of elements such as Y, P, as well as Ca and Mg, which, together with the pH index, suggest an alkaline environment enriched in residues from artificial fertilizers frequently used on arable PAN surfaces. At the other end of the scale, soils formed on clastic rocks (GEO3 in particular) are acidic, revealing a similar mechanism of soil acidification as in the mountainous areas of the DIN region investigated earlier (HASAN et al., 2020). This process, occurring simultaneously with illuviation, may dominate soil formation over the siliciclastic bedrock, particularly in PAN mountainous areas. In contradistinction to DF1, the second discriminant function (DF2), explaining a further 17.97% of the total variability of the Pannonian GEO model (over 50% taken together), indicates a direct relationship between the soils formed on Quaternary deposits and those developed on siliciclastic parent rocks outcropping on the PAN mountains. Grouping of alluvial sediments (GEO13) on the lower part of the diagram, together with metamorphic rocks (GEO1) and mainly basic rocks (GEO2), demonstrates that weathering Table 6. Tests of residual roots (discriminant functions) for all three models. DF Eigen value Eigen (%) Eigen cum Canon. R Wilks’ l c2 df p-level GEOLOGY 1 0.849 33.31 33.31 0.678 0.119 2623.7 336 0.000 2 0.458 17.97 51.28 0.560 0.220 1865.9 297 0.000 3 0.330 12.95 64.23 0.498 0.321 1401.2 260 0.000 4 0.233 9.14 73.37 0.434 0.427 1049.9 225 0.000 5 0.176 6.90 80.27 0.387 0.526 792.1 192 0.000 6 0.149 5.84 86.11 0.360 0.619 592.1 161 0.000 7 0.117 4.59 90.70 0.324 0.711 420.9 132 0.000 8 0.094 3.69 94.39 0.293 0.794 284.2 105 0.000 9 0.057 2.24 96.63 0.232 0.869 173.2 80 0.000 10 0.045 1.76 98.39 0.209 0.918 104.9 57 0.000 11 0.025 0.98 99.37 0.155 0.960 50.1 36 0.060 12 0.016 0.63 100.00 0.127 0.984 20.0 17 0.272 SOIL 1 1.041 51.59 51.59 0.714 0.205 1958.0 196 0.000 2 0.408 20.22 71.81 0.538 0.418 1076.6 162 0.000 3 0.206 10.21 82.02 0.413 0.589 654.1 130 0.000 4 0.171 8.47 90.49 0.382 0.710 423.3 100 0.000 5 0.097 4.81 95.30 0.297 0.831 228.5 72 0.000 6 0.067 3.32 98.62 0.251 0.912 114.4 46 0.000 7 0.028 1.38 100.00 0.165 0.973 34.0 22 0.049 CLC 1 0.322 59.85 59.85 0.494 0.617 597.0 84 0.000 2 0.163 30.30 98.15 0.375 0.816 251.4 54 0.000 3 0.053 9.85 100.00 0.225 0.950 64.1 26 0.000 Table 7. Multivariate test for overall significance of discrimination. Models GEOLOGY SOIL CLC No. of groups 13 8 4 Wilks’ lambda 0.119 0.205 0.617 Approximate F ratio 8.505 11.118 7.628 Degrees of freedom [336; 13463] [196; 8314] [84; 3659] p level p < 0.000 p < 0.000 p < 0.000 G eologia C roatica 165Pjanić et al.: Geochemical Map of Croatia – Key findings and legacy of geochemical mapping in the Pannonian Basin Figure 3. Comparison between variables and groups (centroids) in the discriminant function models GEOLOGY, SOIL and CLC of the clr-transformed data. In the GEOLOGY model: scatterplots of a) variable loadings and b) groups in the reduced discriminant space of the first two discriminant functions (DF1–DF2); in the SOIL model: scatterplots of c) variable loadings and d) DF1–DF2; in the CLC model: scatterplots of e) variable loadings and f) DF1–DF2. G eo lo gi a C ro at ic a 166 Geologia Croatica 78/2 products of the latter are transported to the numerous streams feeding the alluvial plains. Concerning the rest of the Quaternary deposits (GEO10–GEO11–GEO12), as opposed to GEO13 on DF2, the pertinent soils have accumulated a host of elements of low mobility such as Ti, Nb, La and Th, most probably as a result of the high stability of their hosting minerals (both oxides and silicates), such as monazite and xenotime. Also, the occurrence of Zr on DF2, compatible with the soils developed on the latter group of Quaternary deposits, indicates the presence of detrital heavy minerals such as zircon, suggesting either their origin from related DF2 lithological facies, (GEO9–GEO8–GEO4 in particular), or from external material of aeolian origin. Soils developed on GEO13/GEO1– GEO2 lithologies appear enriched mostly in PHE (Potentially Harmful Elements) (Zn, Pb, Hg), possibly of anthropogenic origin. However, the absence of elements constituting the clay minerals (K, Na and Al, together with Ti, Ba and others) calls for attention, and a solution must be sought in residual variance that may explain such unusual behaviour. In effect, DF3 as a true monopolar function (explaining further 12.95% of the total variability) clarifies the unique position of the metamorphic rocks (GEO1) in the overall scheme betraying the latter as a main source of clay­forming elements in the pertinent soils and their acidulous character (low pH; Fig. 3a, b). 3.1.2. SOIL model Soil as the final product of the “state function” inherently includes all the so­called state factors that more or less participated in its formation. Thus, the picture displayed by the SOIL­model DFA scatterplots is the natural canvas painted on the Earth's surface in its most recent geological history. These plots (Fig. 3c, d) represent a collection of soil groups with their associated geochemical imprints originating from various soil­forming processes characteristic of the Croatian PAN region. DF1 explains almost 52% of the total variability of the soil system emphasizing the ubiquity of Fluvisols (FL), the soil group characteristic for the investigated region, since its formation is related to environments that are, by definition, periodically flooded areas of alluvial plains, river fans and valleys (Fig. 2c, d; Fig. 4c, d; Tables 2, 6). Naturally, its formation is fixed on recent (mainly Holocene) f luvial, lacustrine or alluvial deposits, the source material of which originates from rocks comprising the inner Pannonian Mountains (as defined in the GEOLOGY model). In this regard, one cannot neglect the alkaline character of PAN Fluvisols with high pH and enrichment in elements such as Ca and Mg. However, at the same time, it is obvious that these soils are depots of various PHEs, including amongst others, Pb, Zn, Hg, Cu, Ni, which enter the system through various human activities, mostly industry (ROMIĆ & ROMIĆ, 2003; ŠORŠA et al., 2018a, b). In addition, the presence of P, although to a lesser degree, indicates agricultural activity on Fluvisols, affecting their capacity to sustain plant or animal productivity in the long term. Relatively close to FL stands the group of urban soils (URB) where the problem with PHEs is further highlighted. At the other end of the spectrum, Podzols (PZ) and Luvisols (LV) are enriched in clayey components and characterised by their excessive acidity. Contrary to DF1, the second discriminant function (DF2), explaining 20.22% of the total variability is more obviously bipolar contrasting Leptosols (LP) and, partly, Cambisols (CM) with Gleys (GL) and Chernozem (CH), based on the enrichment/depletion relationship of P and Y (present in phosphogypsum fertilizers) on the one side, and Mn on the other. This picture separates agricultural from other soils, since yttrium is widely applied in modern industry and farming, while phosphorus is one of the essential nutrients, occurring in agricultural areas used for permanent crops (phosphogypsum fertilizers). However, phosphorus is one of the main factors of eutrophication, especially in low­lying areas with poor drainage, where dissolved nutrients may stimulate the growth of aquatic species, resulting in anoxia. This process is particularly relevant to Gleys (GL) saturated with (oxygen­depleted) groundwater. In contrast, Leptosols (LP) are generally well­drained soils, not suitable for agricultural utilisation but with characteristics favourable for forestry. In this regard, the role of Mn as a critical micronutrient for plants may be accumulated in surface soils (Leptosols in this case) after a few decades as a consequence of cycling through vegetation, that is, uptake and litterfall. As for the third discriminant function (DF3, 10.21%), it is also concerned with Leptosols but in this case focuses on the distribution of heavy metals (Cu, Ni, Co) that are most probably related to the underlying ophiolitic bedrock building the central PAN area. 3.1.3. CLC model The most typical level (Label 1) of the CLC model, generally relating to the vegetation and animal properties and ecosystem as a whole in Jenny’s extended soil functional­factorial model. JENNY (1941), as described in BOCKHEIM & GENNA­ DIYEV (2010), explains over 90% of the total variability by the first two (DF1 and DF2) of the three discriminant functions altogether (Fig. 3e, f; Fig. 4e, f; Table 6). The first discriminant function is essential (59.85%). It highlights the strong contrast between soils developed on agricultural and artificial surfaces (AGRS/ARTS) with those formed in forests and semi­natural areas and wetlands, when taken together (FSNA/WETL). This distinction is based essentially on the strong presence of the same elements as in the GEO model, referring to groups of Quaternary lithology (GEO10–GEO13), namely P, Y (Ca) and high pH coefficient revealing the alkaline conditions. A similar impact of the inescapable anthropogenic impact of fertilizers is expected between these GEO and CLC groups. However, in the latter case, Cu also added to AGRS an ARTS group (Fig. 3e, f), indicating that a portion of arable land may be covered with vineyards (cupric sulphate) (ROMIĆ et al., 2004). The second CLC discriminant function (DF2) (additional 30.30% of total variance) is strongly monopolar, separating wetland surfaces from all other groups. The WETL group is generally alkaline, characterised by the presence of Ca and Mg, including a high pH index (most probably from sewage systems). However, the occurrence of PHE such as Zn, Hg and Pb indicates that lowlands are highly susceptible to heavy metal accumulation additionally affecting the ecosystem by increasing toxicity, especially along the regional rivers such G eologia C roatica 167Pjanić et al.: Geochemical Map of Croatia – Key findings and legacy of geochemical mapping in the Pannonian Basin Figure 4. Discriminant score maps representing areal distribution of the first and second discriminant function of the models: GEOLOGY – a) DF1, b) DF2; SOIL – c) DF1, d) DF2; and CLC – e) DF1, f) DF2. The increasing influence of the respective geochemical signatures is displayed in warm colours (yellow-orange-red). G eo lo gi a C ro at ic a 168 Geologia Croatica 78/2 as the Drava and Sava (PAVLOVIĆ et al., 2004; PAVLOVIĆ et al., 2019; OREŠČANIN et al., 2004; ŠAJN et al., 2011; HALAMIĆ et al., 2012). 3.1.4. Classification issues Classification efficacy serves as a robust indicator by which the integrity of the previously (a priori) delineated groups can be checked, weighing mathematically “predicted” (computed) against original (observed) classifications (Table 5). This reasoned evaluation of group membership of each sample reveals that a considerable proportion of the original affiliations proved incorrect, which is, without doubt, a result of the complex interaction between state factors involved in soil formation. In contradistinction to the DIN region (HASAN et al., 2020), where the preponderance of the carbonate bedrock strongly influenced the active soil forming processes toward the characteristic patterns in each investigated model, resulting in over 70% correct assignments in each division, the PAN region is obviously more complex. Unquestionably, the increased variation in bedrock lithologies can be considered the main cause of this image as the diverse siliciclastic bedrock, composed of a variety of sedimentary, igneous and metamorphic rocks, is not easy to separate based solely on geochemistry. It is particularly typical for the youngest, Quaternary sediments, where, for example, the majority of soil samples collected over the units mapped as deluvial­proluvial deposits (GEO12) changed their group membership to loess and marsh deposits (GEO10), resulting in a mere 5.26% correct affiliation (Table 8). Incidentally, it is a typical gravity­centring pattern of group assignment where unstable groups tend to centre their “lost” members on some group with the highest number of samples. The analogue in the Dinaric (DIN) region is, quite naturally, the carbonate lithology which is why most of the samples that were collected on carbonate bedrock were labelled á priori as “carbonate clastic rocks” assembled post hoc around the more general lithological group defined as “carbonate rocks in general” (HASAN et al., 2020). In the SOIL model, the samples classified as urban soils (URB), with merely 7.69% of correct memberships, changed Table 8. Classification matrix. Observed groups Predicted groups GEOLOGY G1 G2 G3 G4 G5 G6 G7 G8 G9 G10 G11 G12 G13 Total Correct % G1 21 0 2 0 1 1 0 0 3 3 0 0 1 32 65.63 G2 0 4 0 0 0 2 1 0 2 0 0 0 0 9 44.44 G3 1 0 16 0 0 1 2 2 0 0 0 0 0 22 72.73 G4 0 0 2 2 1 0 0 0 0 2 0 0 0 7 28.57 G5 1 0 1 0 7 0 1 0 1 3 0 0 2 16 43.75 G6 2 0 0 0 0 6 0 0 2 0 0 0 2 12 50.00 G7 1 0 1 0 0 2 12 5 3 5 0 1 2 32 37.50 G8 1 0 2 1 0 2 1 30 18 18 1 0 3 77 38.96 G9 3 0 7 0 1 1 1 7 73 46 2 3 12 156 46.79 G10 5 0 0 0 0 1 2 4 31 394 6 2 44 489 80.57 G11 0 0 0 0 0 0 0 0 0 14 13 1 4 32 40.63 G12 1 0 1 0 0 0 0 2 5 25 0 2 2 38 5.26 G13 0 2 3 0 3 1 1 7 6 101 10 1 197 332 59.34 Total 36 6 35 3 13 17 21 57 144 611 32 10 269 1254 61.96 SOIL FL CM PZ GL LV LP URB CH Total Correct % FL 104 1 0 33 2 3 1 1 145 71.72 CM 12 74 21 18 34 8 2 13 182 40.66 PZ 0 20 104 54 66 4 2 0 250 41.60 GL 18 7 39 178 34 6 2 9 293 60.75 LV 3 23 49 32 138 7 0 5 257 53.70 LP 2 19 4 5 20 33 0 0 83 39.76 URB 3 1 1 3 3 1 1 0 13 7.69 CH 1 0 0 2 3 0 0 25 31 80.65 Total 143 145 218 325 300 62 8 53 1254 52.39 CLC FSNA AGRS ARTS WETL Total Correct % FSNA 284 186 6 29 505 56.24 AGRS 107 536 7 12 662 80.97 ARTS 2 39 7 1 49 14.29 WETL 1 19 0 18 38 47.37 Total 394 780 20 60 1254 67.38 G eologia C roatica 169Pjanić et al.: Geochemical Map of Croatia – Key findings and legacy of geochemical mapping in the Pannonian Basin their affiliation by being deployed to various other soil types, mostly Gleysols (GL), Luvisols (LV) and Fluvisols (FL) indicating, perhaps, inappropriate labelling for this type of soil which denotes more environmental (human­induced, in particular) than soil­based characteristics. This implies that the unmapped soils in the urban area are not predominantly under anthropogenic influence, but retain their pedological properties. As the number of sampling sites in the URB group is small (13 sites), this conclusion is rather tentative. In bypassing the natural soil­forming processes in this way, this label becomes devoid of true meaning in the soil classification system (Table 8). In the CLC model with only four á priori defined groups, the misclassification is emphasised further, especially with regards to the group of artificial surfaces (ARTS), the classification rate of which amounts to a mediocre 14.29%. Similarly, as in the case of “urban soils” (URB) in the SOIL model that lose their members to more “natural” groups, here the ARTS division distributes its members to the group of agricultural surfaces (AGRS) (Tables 5, 8) indicating instability of the so­called “human­impacted” classes whenever they are used for classification purposes. As for the wetland surfaces (WETL), half of its members change their affiliation to agricultural surfaces (AGRS) emphasizing the ubiquitous importance of this latter type of environment in the Pannonian region, specifically in the areas covered by overbank and floodplain sediments. 4. CONCLUSIONS The application of DFA to geological, soil, land use and land cover data from the Pannonian Basin provides significant insights into the geochemical and environmental processes in the topsoil layer. The results highlight lithological and pedological diversity, and consequently land use and land cover, with the first two discriminant functions DF1 and DF2 explaining a significant proportion of the total variance. 1. In this study, discriminant function models were built based on compositional data analysis for the geochemical survey of topsoils in the Pannonian Basin of Croatia and cross­comparison of their geochemical signatures among the selected groups of geological, pedological and land use/land cover divisions. This approach, combining DFA­derived geochemical models with comprehensive environmental data, provides clear insights into the relationships between topsoil processes and geochemical variations according to geology, pedology, land use and land cover. 2. The GEOLOGY model with the first discriminant function – DF1 explains one third of the variance (33.31%) and shows the clear bipolar relationship between Quaternary sediments and older rocks. This function distinguishes topsoils on Quaternary deposits, which are enriched in elements indicative of an alkaline environment from soils formed on older, more acidic rocks such as ophiolites and clastic rocks containing some carbonates. The second discriminant function – DF2 clarifies 17.97% of the total variability. DF2 also discriminates between soils on alluvial sediments and those on metamorphic and ophiolitic rocks, indicating a significant transport of weathering products in alluvial plains. DF3 explains 12.95% of the variability and emphasises the unique geochemical signature of soils on metamorphic rocks and their role in the formation of clay minerals and soil acidification. 3. The SOIL model reflects the multi­layered nature of soil formation processes, with DF1 accounting for the largest proportion of variability (51.59%). This feature emphasises the prevalence of Fluvisols in the region, associated with recent fluvial and alluvial deposits, as well as their enrichment with essential nutrients and potentially harmful elements due to human activities. DF2 explains 20.22% of the total variability and is almost bipolar. It compares agricultural soils enriched with phosphorus and yttrium from fertilisers with other soil types and shows the impact of agricultural practices on soil chemistry. DF3 describes only 10.21% of the variability and focuses on the distribution of heavy metals in Leptosols, especially those overlying ophiolitic bedrock, and explains the influence of underlying geology on soil composition. 4. The CLC model with DF1 (59.85%) and DF2 (30.30%) together explaining over 90% of the total variability, describes the strong contrast between soils on agricultural/artificial land and those in forests, semi­ natural areas and wetlands. DF1 highlights the anthropogenic inf luences on agricultural soils characterised by the presence of phosphorus, yttrium and high pH, while DF2 emphasizes wetlands characterised by alkaline conditions and the accumulation of heavy metals, ref lecting their susceptibility to pollution and the accumulation of heavy metals from regional rivers. 5. The effectiveness of the classification shows that the complex factors of soil formation led to considerable misclassification. The diversity of lithological units, especially in Quaternary sediments, makes it difficult to delineate geochemical groups, for example, deluvial­ proluvial deposits (GEO12) only had an affiliation rate of 5.26%. Urban soils, with only 7.69% being a correct classification, are often misclassified as other soil types, suggesting that their current classification does not accurately reflect their pedological properties. In addition, artificial surfaces had a high misclassification rate of 14.29%. They were frequently assigned to agricultural land, indicating anthropogenic influence. The results of this study illustrate the geochemical diversity and the complex interactions between geology, soil formation processes, land use and land cover in the Pannonian Basin. The geochemical signatures of geology, lithology and land use provide an authentic framework for environmental management and further research. ACKNOWLEDGEMENT This study was financially supported by The Ministry of Science, Education and Sports, Republic of Croatia – Scientific Project: The Basic Geochemical Map of Croatia, Number 181­ G eo lo gi a C ro at ic a 170 Geologia Croatica 78/2 1811096­1181. We are grateful to the Associate Editor academic Goran Durn and the anonymous reviewers for their thorough evaluations which significantly improved the manuscript. REFERENCES ACME Analytical Laboratories Ltd. (2007): Assaying and geochemical analy­ ses.– ACME Analytical Laboratories Ltd., Vancouver, 19 p. AITCHISON, J. (1986): The Statistical Analysis of Compositional Data.– Cha­ pman and Hall, London, New York, 416 p. http://doi.org/10.1007/978­94­ 009­4109­0 AITCHISON, J. (2008): The single principle of compositional data analysis, continuing fallacies, confusions and misunderstandings and some sugge­ sted remedies.– In: DAUNIS­I­ESTADELLA, J. & MARTÍN­FERNÁN­ DEZ, J. (eds.): Proceedings of CoDaWork’08, The 3rd Compositional Da­ ta Analysis Workshop. University of Girona, Girona, 1–27. AITCHISON, J. & EGOZCUE, J.J. (2005): Compositional data analysis: whe­ re are we and where should we be heading?– Mathematical Geology, 37/7, 829–850. http://doi.org/10.1007/s11004­005­7383­7 BARCELÓ­VIDAL, C., MARTIN­FERNANDEZ, J.A. & PAWLOWSKY­ GLAHN, V. (1999): Comment on „Singularity and Nonnormality in the Classification of Compositional Data“ by G.C. Bohling, J.C. Davis, R.A. Olea & J. Harff.– Mathematical Geology, 31, 581–585. BAŠIĆ, F. (2013): The Soils of Croatia.– Springer Dordrecht, 179 p. https://doi. org/10.1007/978­94­007­5815­5 BAŠIĆ, F., BOGUNOVIĆ, M., BOŽIĆ, M., HUSNJAK, S., JURIĆ, I., KISIĆ, I., MESIĆ, M., MIROŠEVIĆ, N., ROMIĆ, D. & ŽUGEC, I. (2007): Re­ gionalisation of Croatian Agriculture.– Agriculturae Conspectus Scien­ tificus, 72/1, 27–38. BOCKHEIM, J.G. & GENNADIYEV, A.N. (2010): Soil­factorial models and earth system science: a review.– Geoderma, 159/3–4, 243–251. https://doi. org/10.1016/j.geoderma.2010.09.005 BOCKHEIM, J.G., GENNADIYEV, A.N., HAMMER, R.D. & TANDARICH, J.P. (2005): Historical development of key concept in pedology.– Geoder­ ma, 124, 23–36. https://doi.org/10.1016/j.geoderma.2004.03.004 BOGUNOVIĆ, M., VIDAČEK, Ž., RACZ, Z., HUSNJAK, S., ŠPOLJAR, A. & SRAKA, M. (1998): FAO UNESCO pedological map of the Republic of Croatia 1:1.000.000.– University of Zagreb, Faculty of Agriculture, Zagreb. BRUNOVIĆ, D., MIKO, S., ILIJANIĆ, N., PEH, Z., HASAN, O., KOLAR, T., ŠPARICA MIKO, M. & RAZUM, I. (2019): Holocene foraminiferal and geochemical records in the coastal karst dolines of Cres Island, Croatia.– Geologia Croatica, 72/1, 19–42. https://doi:10.4154/gc.2019.02 BUCCIANTI, A. (2013): Is compositional data analysis a way to see beyond the illusion?– Computers & Geosciences, 50, 165–173. http://doi.org/10. 1016/j.cageo.2012.06.012 BUCCIANTI, A. & GRUNSKY, E. (2014): Compositional data analysis in geo­ chemistry: are we sure to see what really occurs during natural proces­ ses?– Journal of Geochemical Exploration, 141, 1–5. https://doi.org/10. 1016/j.gexplo.2014.03.022 DAUNIS­I­ESTADELLA, J., THIÓ­HENESTROSA, S. & MATEU­FIGUE­ RAS, G. (2011): Two more things about compositional biplots: quality of projection and inclusion of supplementary elements.– In: EGOZCUE, J.J., TOLOSANA­DELGADO, R. & ORTEGO, M.I. (eds.): Proceedings of the 4th International Workshop on Compositional Data Analysis (Co­ DaWork'11). University of Girona, Girona, 1–14. DAVIS, J.C. (1973): Statistics and Data Analysis in Geology.– John Wiley & Sons, New York. 550 p. DAVIS, J.C. (1986): Statistics and Data Analysis in Geology, 2nd Edition.– Jo­ hn Wiley & Sons, New York, 656 p. DILLON, W.R. & GOLDSTEIN, M. (1984): Multivariate analysis: Methods and applications.– John Wiley & Sons, New York, 608 p. EGOZCUE, J.J. & PAWLOWSKY­GLAHN, V. (2006): Simplicial geometry for compositional data.– In: BUCCIANTI A., MATEU­FIGUERAS, G. & PAWLOWSKY­GLAHN, V. (eds.): Compositional Data Analysis in the Geosciences: From Theory to Practice. Geological Society London, Special Publication, 264, 145–158. http://doi.org/10.1144/GSL.SP.2006. 264. 01.11 FILIPČIĆ, A. (2023): Klima Hrvatske.– In: MAGAŠ, D (ed.): Velika geografi­ ja Hrvatske, knjiga 2. – Fizička geografija Hrvatske – Prirodno­geograf­ ska osnova razvoja [Major Geography of Croatia, Book 2 – Psyhical Ge- ography of Croatia – in Croatian], University of Zadar, Department of Geography, Školska knjiga, Zagreb, Zadar, 432–467 p. GALOVIĆ, L. & PEH, Z. (2016): Mineralogical discrimination of the Pleisto­ cene loess/paleosol sections in Srijem and Baranja, Croatia.– Aeolian Re­ search, 21, 151–162. http://doi.org/10.1016/j.aeolia.2016.04.006 GALOVIĆ, L., HUSNJAK, S., ŠORŠA, A. & MARTINČEVIĆ LAZAR, J. (2023): Evidence and mineralogical and physico­chemical properties of chernozem and chernozem­like soils in Croatia.– Geologia Croatica, 76/3, 113–129. https://doi.org/10.4154/gc.2023.08 GALOVIĆ, L., ŠORŠA, A. & PEH, Z. (2024): Modelling the paleopedological development of loess/palaeosol sections using discriminant function analysis of compositional geoscience data, case study in Eastern Croatia.– Journal of Geochemical Exploration, 262, 107473. https://doi.org/10.1016/j. gexplo.2024.107473 GIZDAVEC, N., GAŠPAROVIĆ, M., MIKO, S., LUŽAR­OBERITER, B., ILI­ JANIĆ, N. & PEH, Z. (2022): Discrimination of Rock Units in Karst Terrains Using Sentinel­2A Imagery.– Remote Sensing, 14/20, 5169. http://doi:10.3390/rs14205169 GRIZELJ, A., PEH, Z., TIBLJAŠ, D., KOVAČIĆ, M. & KUREČIĆ, T. (2017): Mineralogical and geochemical characteristics of Miocene pelitic sedi­ mentary rocks from the south­western part of the Pannonian Basin System (Croatia): Implications for provenance studies.– Geoscience Frontiers, 8/1, 65–80. https://doi.org/10.1016/j.gsf.2015.11.009 HALAMIĆ, J. & MIKO, S. (eds.) (2009): Geochemical atlas of the Republic of Croatia.– Croatian Geological Survey, Zagreb, 87 p. HALAMIĆ, J., MIKO, S., PEH, Z. & GALOVIĆ, L. (2000): Upute za izradu "Osnovne geokemijske karte Republike Hrvatske" [Guidelines for geoc- hemical mapping of Croatia – in Croatian]. Fond struč. dok., Croatian Geological Survey, 51/2000, 22 p. HALAMIĆ, J., PEH, Z., MIKO, S., GALOVIĆ, L. & ŠORŠA, A. (2012): Geo­ chemical Atlas of Croatia: environmental implications and geodynamical thread.– Journal of Geochemical Exploration, 115, 36–46. https://doi.or­ g/10.1016/j.gexpl o.2012.02.006 HASAN, O., MIKO, S., ILIJANIĆ, N., BRUNOVIĆ, D., DEDIĆ, Ž., ŠPARICA MIKO, M. & PEH, Z. (2020): Discrimination of topsoil environments in a karst landscape: An outcome of a geochemical mapping campaign.– Geochemical Transactions, 21, 1. https://doi.org/10.1186/s12932­019­ 0065­z HGI­CGS (HRVATSKI GEOLOŠKI INSTITUT­CROATIAN GEOLOGICAL SURVEY) (1992): Geološke karte sjeverozapadne Hrvatske, M 1:100 000, Vrh temeljnog gorja [Geological Maps of Northwestern Croatia, Scale 1:100,000, The Top of the Mesozoic Basement].– Archives of HGI­CGS, 104/92, Zagreb. HGI­CGS (HRVATSKI GEOLOŠKI INSTITUT­CROATIAN GEOLOGICAL SURVEY) (2009): Geološka karta Republike Hrvatske M 1:300 000 [Ge- ological Map of Republic of Croatia, M 1:300,000].– HGI­CGS, Depar­ tment for Geology, Zagreb. https://www.hgi­cgs.hr/geoloske­karte/. ILIJANIĆ, N., KOVAČEVIĆ GALOVIĆ, E., GIZDAVEC, N., IVKIĆ FILIPO­ VIĆ, I., MIKO, S. & PEH, Z. (2023): Geochemical records in subaerial exposure environments in Croatia using discriminant function analysis of bauxite data.– Frontiers in Earth Science, 10, 1055435. https://doi. org/10.3389/feart.2022.1055435 ISO (2002a): ISO 10381­1: Soil quality–sampling–part 1: guidance on the de­ sign of sampling programs.– International Organization for Standardiza­ tion. ISO (2002b): ISO 10381­2: Soil quality–sampling–part 2: guidance on the de­ sign of sampling programs.– International Organization for Standardiza­ tion. JAMIČIĆ, D. & CRNKO, J. (2009): Kompleks metamorfnih stijena Psunja, Krndije i Moslavačke gore, Prekambrij, Pk [Complex of metamorphic rocks of the Psunj, Krndija and Moslavačka Gora Mts., Precambrian, http://doi.org/10.1007/978-94-009-4109-0 http://doi.org/10.1007/978-94-009-4109-0 http://doi.org/10.1007/s11004-005-7383-7 https://doi.org/10.1007/978-94-007-5815-5 https://doi.org/10.1007/978-94-007-5815-5 https://doi.org/10.1016/j.geoderma.2010.09.005 https://doi.org/10.1016/j.geoderma.2010.09.005 https://doi.org/10.1016/j.geoderma.2004.03.004 https://doi:10.4154/gc.2019.02 http://doi.org/10.1016/j.cageo.2012.06.012 http://doi.org/10.1016/j.cageo.2012.06.012 https://doi.org/10.1016/j.gexplo.2014.03.022 https://doi.org/10.1016/j.gexplo.2014.03.022 http://doi.org/10.1144/GSL.SP.2006.264.01.11 http://doi.org/10.1144/GSL.SP.2006.264.01.11 http://doi.org/10.1016/j.aeolia.2016.04.006 https://doi.org/10.4154/gc.2023.08 https://doi.org/10.1016/j.gexplo.2024.107473 https://doi.org/10.1016/j.gexplo.2024.107473 http://doi:10.3390/rs14205169 https://doi.org/10.1016/j.gsf.2015.11.009 https://doi.org/10.1016/j.gexpl o.2012.02.006 https://doi.org/10.1016/j.gexpl o.2012.02.006 https://doi.org/10.1186/s12932-019-0065-z https://doi.org/10.1186/s12932-019-0065-z https://www.hgi-cgs.hr/geoloske-karte/ https://doi.org/10.3389/feart.2022.1055435 https://doi.org/10.3389/feart.2022.1055435 G eologia C roatica 171Pjanić et al.: Geochemical Map of Croatia – Key findings and legacy of geochemical mapping in the Pannonian Basin Pk – in Croatian].– In: VELIĆ, I. & VLAHOVIĆ, I. (eds.): Tumač za Ge­ ološku kartu Hrvatske M 1:300.000 [Explanatory Notes for the Geologi- cal Map of Croatia M 1:300,000– in Croatian].– Croatian Geological Sur­ vey, 11–13. JENNY, H. (1941): Factors of Soil Formation: A System of Quantitative Pedo­ logy.– Dover Publications, New York, 271 p. KÖPPEN, W. (1918): Klassifikation der Klimate nach Temperatur, Niederschlag und Jahreslauf, Petermanns Geographische Mitteilungen, 64, 193–203. KOVAČEVIĆ GALOVIĆ, E., ILIJANIĆ, N., PEH, Z., MIKO, S. & HASAN, O. (2012): Geochemical discrimination of Early Palaeogene bauxites in Croatia.– Geologia Croatica, 65, 53–65. https://doi.org/10.4154/gc.2012.04 MERT, M.C., FILZMOSER, P., ENDEL, G. & WILBACHER, I. (2018): Com­ positional data analysis in epidemiology.– Statistical Methods in Medical Research, 27/6, 1878–1891. https://doi.org/10.1177/0962280216671536 MIKO, S., HALAMIĆ, J., PEH, Z. & GALOVIĆ, L. (2001): Geochemical ba­ seline mapping of soils developed on diverse bedrock from two regions in Croatia.– Geologia Croatica, 54, 53–118. MUTIĆ, R. (1975): Sedimentološka ispitivanja naslaga lesa iz okolice Vinko­ vaca, Našica i Valpova [Sedimentological studies of loess deposits from the vicinity of Vinkovci, Našice and Valpovo– in Croatian].– Geološki Vjesnik, 28, 269–286. OREŠČANIN, V., LULIĆ, S., PAVLOVIĆ, G. & MIKELIĆ, L. (2004): Granu­ lometric and chemical composition of the Sava River sediments upstream and downstream of the Krško nuclear power plant.– Environmental geo­ logy (Berlin), 46, 605–613. https://doi.org/10.1007/s00254­004­1066­4 PAMIĆ, J. & LANPHERE, M. (1991): Hercynian granites and metamorphic rocks from the Papuk, Psunj, Krndija and the surrounding basement of the Pannonian Basin (Northern Croatia, Yugoslavia).– Geologija, 34, 81– 253. PANAGOS, P. (2006): The European soil database.– GEO: connexion, 5/7, 32–33. PAVELIĆ, D. & KOVAČIĆ, M. (2018): Sedimentology and stratigraphy of the Neogene rift­type North Croatian Basin (Pannonian Basin System, Cro­ atia): A review.– Marine and Petroleum Geology, 91, 455–469. https://doi. org/10.1016/j.marpetgeo.2018.01.026 PAVLOVIĆ, G., PROHIĆ, E. & TIBLJAŠ, D. (2004): Statistical assessment of geochemical pattern in overbank sediments of the river Sava, Croatia.– Environmental Geology, 46, 132–143. https://doi.org/10.1007/s00254­004­ 1019­y PAVLOVIĆ, P., MARKOVIĆ, M., KOSTIĆ, O., SAKANC, S., ĐORĐEVIĆ, D., PEROVIĆ, V., PAVLOVIĆ, D., PAVLOVIĆ, M., ČAKMAK, D., JA­ RIĆ, S., PAUNOVIĆ, M. & MITROVIĆ, M. (2019): Evaluation of poten­ tially toxic element contamination in the riparian zone of the River Sava.– Catena, 174, 399–412. https://doi.org/10.1016/j.catena.2018.11.034 PAWLOWSKY­GLAHN, V. & EGOZCUE, J.J. (2006): Compositional data and their analysis: an introduction.– In: BUCCIANTI, A., MATEU­FI­ GUERAS, G. & PAWLOWSKY­GLAHN. V. (eds.): Compositional Data Analysis in the Geosciences: From Theory to Practice. Geological Society London, Special Publication, 264, 1–10. PAWLOWSKY­GLAHN, V., EGOZCUE, J.J. & TOLOSANA­DELGADO, R. (2007): Lecture notes on compositional data analysis.– Universitat de Gi­ rona, Girona, 96 p. PEH, Z. & HALAMIĆ, J. (2010): Discriminant function model as a tool for classification of stratigraphically undefined radiolarian cherts in ophioli­ te zones.– Journal of Geochemical Exploration, 107, 30–38. https://doi. org/10.1016/j.gexplo.2010.06.003 PEH, Z. & KOVAČEVIĆ GALOVIĆ, E. (2014): Geochemistry of Istrian Lower Palaeogene bauxites—Is it relevant to the extent of subaerial exposure during Cretaceous times?– Ore Geology Reviews, 63, 296–306. https:// doi.org/10.1016/j.oregeorev.2014.05.020 PEH, Z. & KOVAČEVIĆ GALOVIĆ, E. (2016): Geochemistry of Lower Pala­ eogene bauxites–unique signature of the tectonostratigraphic evolution of a part of the Croatian Karst.– Geologia Croatica, 69/2, 269–279. https:// doi.org/10.4154/gc.2016.24 PEH, Z., ŠAJN, R., HALAMIĆ, J. & GALOVIĆ, L. (2008): Multiple discrimi­ nant analysis of the Drava River alluvial plain sediments.– Environmen­ tal Geology, 55/7, 1519–1535. https://doi.org/10.1007/s0025 4­007­1102­2 PEH, Z., MIKO, S. & HASAN, O. (2010): Geochemical background in soils: a linear process domain? An example from Istria (Croatia).– Environmen­ tal Earth Sciences, 59/6, 1367–1383. https://doi.org/10.1007/s12665­009­ 0125­2 PHILLIPS, J.D. (1998): On the relation between complex systems and the fa­ ctorial model of soil formation.– Geoderma, 86, 1–21. https://doi. org/10.1016/S0016­7061(98)00054­8 PHILLIPS, J.D. (2002): Global and local factors in earth surface systems.– Eco­ logical Modelling, 149/3, 257–272. https://doi.org/10.1016/S0304­ 3800(01)00461­6 PHILLIPS, J.D. (2017): Soil Complexity and Pedogenesis.– Soil Science, 182, 117–127. https://doi.org/10.1097/SS.0000000000000204 POCH, R.M., GALOVIĆ, L., HUSNJAK, S., MARTINČEVIĆ LAZAR, J., HEĆEJ, N., RUŽIČIĆ, S., PJANIĆ, A., ÁLVAREZ, D., BEERTEN, K., GAJIĆ, R., STEJIĆ, P. & PANDUROV, M. (2024): Soil formation and en­ vironmental reconstruction of a loess­paleosol sequence in Zmajevac, Cro­ atia.– Catena, 247, 108507. https://doi.org/10.1016/j.catena.2024.108507 REIMANN, C., FILZMOSER, P., GARRETT, R.G. & DUTTER, R. (2008): Statistical data analysis explained: Applied environmental statistics with R.– John Wiley & Sons, New York, 368 p. http://doi.org/10.1002/ 9780470987605 REIMANN, C., BIRKE, M., DEMETRIADES, A., FILZMOSER, P. & O’CONNOR, P. (eds.) (2014): Chemistry of Europe’s agricultural soils – Part A. Methodology and interpretation of the GEMAS data set.– Geolo­ gisches Jahrbuch (Reihe B 102), Schweizerbarth, Hannover, 528 p. ROCK, N.M.S. (1988): Lecture Notes in Earth Sciences, Vol 18.– In: Numeri­ cal Geology. Springer Berlin, Heidelberg, 427 p. https://doi.org/10.1007/ BFb0045143 ROMIĆ, M. & ROMIĆ, D. (2003): Heavy metals distribution in agricultural topsoils in urban area.– Environmental Geology, 43, 795­805. https://doi. org/10.1007/s00254­002­0694­9 ROMIĆ, M., ROMIĆ, D., DOLANJSKI, D. & STRIČEVIĆ, I. (2004): Heavy Metals Accumulation in Topsoils from the Wine­Growing Regions, Part 1. Factors which Control Retention.– Agriculturae Conspectus Scientifi­ cus, 69/1, 1–10. STARIJAŠ, B., GERDES, A., BALEN, D., TIBLJAŠ, D. & FINGER, F. (2010): The Moslavačka Gora crystalline massif in Croatia: a Cretaceous heat dome within remnant Ordovician granitoid crust.– Swiss Journal of Ge­ osciences, 103/1, 61–82. https://doi.org/10.1007/s00015­010­0007­3 STATSOFT INC. (2006): STATISTICA (Data Analysis Software System), Version 7.1. www.statsoft.com. STRAHLER, A.N. (1980): Systems Theory in Physical Geography.– Physical Geography, 1/1, 1–27. https://doi.org/10.1080/02723646.1980.10642186 ŠAJN, R., HALAMIĆ, J., PEH, Z., GALOVIĆ, L. & ALIJAGIĆ, J. (2011): As­ sessment of the natural and anthropogenic sources of chemical elements in alluvial soils from the Drava River using multivariate statistical met­ hods.– Journal of Geochemical Exploration, 110/3, 278­289, https://doi. org/10.1016/j.gexplo.2011.06.009 ŠORŠA, A. (2009): GIS – data processing and graphical presentation.– In: HALAMIĆ, J. & MIKO, S. (eds.): Geochemical Atlas of the Republic of Croatia.– Croatian Geological Survey, Zagreb, 29. ŠORŠA, A., PEH, Z. & HALAMIĆ, J. (2018a): Geochemical mapping the ur­ ban and industrial legacy of Sisak, Croatia, using discriminant function analysis of topsoil chemical data.– Journal of Geochemical Exploration, 187, 155–167. https://doi.org/10.1016/j.gexplo.2017.07.014 ŠORŠA, A., MILER, M., GOSAR, M. & HALAMIĆ, J. (2018b): Follow­up geochemical studies and mineralogical investigations by scanning ele­ ctron microscopy/energy dispersive spectroscopy (SEM/EDS) of soil samples from the industrial zone of Sisak, Croatia.– Journal of Geoche­ mical Exploration, 187, 168–183. https://doi.org/10.1016/j.gexplo.2017. 07.015 https://doi.org/10.4154/gc.2012.04 https://doi.org/10.1177/0962280216671536 https://doi.org/10.1016/j.marpetgeo.2018.01.026 https://doi.org/10.1016/j.marpetgeo.2018.01.026 https://doi.org/10.1007/s00254-004-1019-y https://doi.org/10.1007/s00254-004-1019-y https://doi.org/10.1016/j.catena.2018.11.034 https://doi.org/10.1016/j.gexplo.2010.06.003 https://doi.org/10.1016/j.gexplo.2010.06.003 https://doi.org/10.1016/j.oregeorev.2014.05.020 https://doi.org/10.1016/j.oregeorev.2014.05.020 https://doi.org/10.4154/gc.2016.24 https://doi.org/10.4154/gc.2016.24 https://doi.org/10.1007/s0025 4-007-1102-2 https://doi.org/10.1007/s12665-009-0125-2 https://doi.org/10.1007/s12665-009-0125-2 https://doi.org/10.1016/S0016-7061(98)00054-8 https://doi.org/10.1016/S0016-7061(98)00054-8 https://doi.org/10.1016/S0304-3800(01)00461-6 https://doi.org/10.1016/S0304-3800(01)00461-6 https://doi.org/10.1097/SS.0000000000000204 https://doi.org/10.1016/j.catena.2024.108507 http://doi.org/10.1002/9780470987605 http://doi.org/10.1002/9780470987605 https://doi.org/10.1007/s00254-002-0694-9 https://doi.org/10.1007/s00254-002-0694-9 https://doi.org/10.1007/s00015-010-0007-3 http://www.statsoft.com https://doi.org/10.1080/02723646.1980.10642186 https://doi.org/10.1016/j.gexplo.2011.06.009 https://doi.org/10.1016/j.gexplo.2011.06.009 https://doi.org/10.1016/j.gexplo.2017.07.014 https://doi.org/10.1016/j.gexplo.2017.07.015 https://doi.org/10.1016/j.gexplo.2017.07.015 G eo lo gi a C ro at ic a 172 Geologia Croatica 78/2 TOLOSANA­DELGADO, R. (2012): Uses and misuses of compositional data in sedimentology.– Sedimentary Geology, 280, 60–79. https://doi.or­ g/10.1016/j.sedgeo.2012.05.005 VELIĆ, I. & VLAHOVIĆ, I. (eds.) (2009): Tumač Geološke karte Republike Hrvatske 1: 300.000 [Explanatory notes for Geological Map of the Republic of Croatia, Scale 1:300,000– in Croatian].– Croatian Geological Survey, 141 p., Zagreb. ZANINOVIĆ, K., GAJIĆ­ČAPKA, M., PERČEC TADIĆ, M., VUČETIĆ, M., MILKOVIĆ, J., BAJIĆ, A., CINDRIĆ, K., CVITAN, L., KATUŠIN, Z., KAUČIĆ, D., LIKSO, T., LONČAR, E., LONČAR, Ž., MIHAJLOVIĆ, D., PANDŽIĆ, K., PATARČIĆ, M., SRNEC, L. & VUČETIĆ, V. (2008): Klimatski atlas Hrvatske 1961–1990, 1971–2000. [Climate atlas of Cro- atia 1961–1990, 1971–2000– in Croatian].– Croatian Meteorological and Hydrological Service, Zagreb, 200 p. Web sources: CORINE LAND COVER (2018): CLC2018 raster data of the Pannonian ba­ sin.– European Environment Agency (EEA). Available at: http://www. eea.europa.eu. Accessed on: April 1st, 2024. https://doi.org/10.1016/j.sedgeo.2012.05.005 https://doi.org/10.1016/j.sedgeo.2012.05.005 http://www.eea.europa.eu http://www.eea.europa.eu