Three-Dimensional Geologic Modeling of the Kiruna Mining District, Sweden: Insights into the Crustal Architecture and Structural Controls on Iron Oxide-Apatite Mineralization

Northern Sweden is a region with significant geologic and economic importance for the European Union and hosts several large mineral deposits, identified and exploited from the 17^th century onward. The Kiruna mining district belongs to the northern Norrbotten ore province and is known as the type locality for Kiruna-type iron oxide-apatite (IOA) deposits, hosting the world-class Kiirunavaara deposit (Geijer, 1919; Parák, 1975; Forsell, 1987). The Kiirunavaara mine is the largest underground iron ore mine in the world, and the ore reserves in the Kiruna area are estimated at 830 million metric tonnes (Mt) with a grade of 41.5% Fe (Luossavaara-Kiirunavaara AB, 2022). Recent exploration efforts indicated resources with a significant amount of apatite-hosted phosphorus and rare earth elements (REEs) in the stratigraphically overlying magnetite-hematite–dominated Per Geijer IOA deposits (Luossavaara-Kiirunavaara AB, 2022). The mining district also hosts other important deposit types, such as the iron oxide copper-gold (IOCG)-style Rakkurijärvi and Pahtohavaare deposits and stratiform copper deposits in the Viscaria and Eastern Pahtohavare area (Martinsson, 1997; Logan et al., 2022).

The Kiruna mining district has been the focus of numerous geologic, geochemical, and geophysical studies during the last century, which resulted in extensive mapping and drilling campaigns. The basin architecture was studied by several authors (e.g., Frietsch, 1979; Martinsson, 2004; Andersson et al., 2021), which resulted in conceptual geologic profiles and models. These models and profiles do not offer a 3-D view of the subsurface structural setting and, due to their conceptual nature, it is hard to further integrate them with newly collected data sets. Geometric 3-D geologic models can help in testing conceptual models against the structural style, geologic framework, and geophysical characteristics based on the data available at a specific temporal instance. Therefore, multiple models can be generated as new data becomes available and tested against up-to-date knowledge (Corbel and Wellmann, 2015). Validated geologic hypotheses based on the 3-D geospatial environment can be used for key economic decisions, such as exploration targeting (de Kemp et al., 2016). However, 3-D camp-scale models must integrate regions of sparse data distributions with regions where abundant data sets (such as local clusters of measurements in densely drilled areas) are available (de Kemp et al., 2016). This approach of camp-scale 3-D modeling has been applied by Schetselaar et al. (2017) in the Flin Flon mining district in Canada and similarly by de Kemp et al. (2016) in regional scale at the Purcell anticlinorium in the East Kootenay region of Canada. In northern Sweden, the crustal architecture and mineralization styles of the Skellefte district (Bauer et al., 2015), Kristineberg area (Skyttä et al., 2013), the Malmberget area (Bauer et al., 2018), and the Nautanen deformation zone (Bauer et al., 2022) have been studied and resulted in 3-D geologic models. Combined geologic and geophysical modeling of the Paleoproterozoic Vakko-Kovo greenstone belts, a few kilometers north of Kiruna, showed the 3-D structure of the crust (Luth et al., 2018). Deposit- to regional-scale 3-D and 4-D models in the Skellefte district, northern Sweden, improved the understanding of the structural framework and geologic history of the district by integrating the deformation history, structural control on sedimentation, volcanism, and mineralization, providing a framework for future industrial and academic activities (Bauer et al., 2015).

Constraining tectonic, metamorphic, and deformation events in Precambrian shield areas with polyphase tectonic histories can be challenging because of the polydeformed state of the bedrock that often experienced additional overprinting thermal events (Royer et al., 2015; Schetselaar et al., 2017; Lahtinen et al., 2018). In northern Sweden, this is further complicated by limited outcrop exposures, bedrock surfaces abraded by repeated glaciations, and flat topography. Even though the Kiruna area is one of the most active mining districts in Sweden, 3-D geologic models from the area focused mainly on local scales, investigating the geometry of the mineralization and to a lesser extent lithological and structural features that affect the mineralization and mining activities. The scientific community has not yet benefited from a 3-D model of the whole mining district that includes the whole succession of supracrustal rocks and in which surface and subsurface observations and measurements are integrated into one model.

In this contribution, we aim to build a 3-D geologic model of the Kiruna mining district to provide insights into geometries and relationships of the structures that are not evident on 2-D maps and cross sections by identifying structural trends and lithological patterns that help in assessing the potential for undiscovered deposits. This study complements and tests hypotheses presented in earlier studies based on surface mapping (e.g., Andersson et al., 2021), investigating the geologic architecture of the Kiruna mining district. The integrated data set includes the digital elevation model, updated geologic map (Fig. 1), aeromagnetic anomaly map (App. Fig. A1), orebody models, drill core logs (including structures, lithology, and geotechnical parameters), outcrop observations, and surface structural measurements. To consider disparities in data resolution and to connect the heterogeneous spatial data, interpretive support is needed; therefore, the elements of the models are variably constrained—from areas constrained only by a few markers or measurements to areas where thousands of markers or measurements are available. Carmichael and Ailleres (2016) highlighted that data collected by geologists can be of large extent and spatially biased, and statistical analysis is needed to determine the trends and to efficiently upscale the data for modeling. Structural analysis was conducted similar to grid cell averaging as an upscaling method during modeling, and structural domains are separated by major structural elements such as faults and shear zones, resulting in a coarser initial distribution of the input data in the Kiruna area. A hybrid 3-D modeling technique is used in this study, which aims at producing a wider spatial perspective of the subsurface. Deposit-scale structural analysis and investigation of geotechnical parameters from drill core measurements focused on identifying the main structures that played an important role during deformation and on the character of the IOA orebodies and their relationship with the hosting horizons.

The oldest rocks in northern Norrbotten are composed of granitic, tonalitic, and amphibolitic gneisses and minor metasedimentary rocks that are part of the Norrbotten craton, which constitutes the basement rocks in this part of the Fennoscandian Shield (2.9–2.6 Ga; Gaál and Gorbatschev, 1987). The southern extent of the basement has been delineated at depth based on Nd isotope signatures in Paleoproterozoic plutonic rocks diagnostic of Archean protoliths (Öhlander et al., 1993; Mellqvist, 1999) and coincides with the Luleå-Jokkmokk zone (illustrated on the inset map in Fig. 1). A period of rifting and associated mafic magmatism between 2.5 and 2.1 Ga resulted in continental breakup of the Karelian craton into three microcontinents including the Norrbotten craton (2.06 Ga; Lahtinen et al., 2005; Skyttä et al., 2019). This continental-scale rifting resulted from E-W– or NE-SW–directed extension, forming a NW-SE–trending Rhyacian tholeiite province stretching from Russia to Norway with associated rift-parallel normal fault systems (Melezhik and Hanski, 2013; Bingen et al., 2015; Skyttä et al., 2019). The importance of rift-related basement structures has been described in northern Finland in the Peräpohja belt interpreted as deep-seated structures controlling the deposition of supracrustal rocks and their subsequent deformation (Piippo et al., 2019; Skyttä et al., 2019). The Archean basement rocks of the Norrbotten craton are overlain by rift-related Rhyacian tholeiitic volcanic and volcaniclastic rocks and sedimentary successions that form ~NNW-SSE– to NNE-NNW–trending greenstone belts in northern Sweden (Pharaoh and Pearce, 1984; Lahtinen et al., 2005; Bingen et al., 2015). The Kovo Group comprising tholeiitic basalt, calc-alkaline volcaniclastic rocks, and clastic sediments forms the base for the Kiruna Greenstone Group and lies unconformably on the Archean rocks (Martinsson, 1997), with a minimum depositional age of 2.18 Ga (Skiöld, 1986). The Kiruna Greenstone Group is composed of a thick pile of mafic to ultramafic rocks (Martinsson et al., 2016). In northern Sweden, the influence of basement structures on the deposition of the Orosirian supracrustal rocks is less understood. However, it has been postulated that preexisting Rhayacian rift structures intersecting the Archean basement played an important role during the Orosirian crustal evolution in northern Norrbotten (Andersson et al., 2021), which is in line with results presented in geophysical and geologic modeling of the Vakko-Kovo greenstone belt north of Kiruna (Luth et al., 2018).

In northern Norrbotten, the change from continental breakup to a period of reassembly of the dispersed cratons involved subduction and accretion toward the northeast along an active continental margin (Bergman et al., 2001; Lahtinen et al., 2005; Skyttä et al., 2019). Remnants of the northeasterly-dipping subducted slab were identified approximately 350 km south of Kiruna by using magnetotellurics (Rasmussen et al., 1987) and seismic reflection profiles (BABEL Working Group, 1990). The transition marked the start of the polyphase Svecokarelian orogenic cycle (1.90–1.78 Ga) and started with the development of an intracontinental back-arc basin (ca. 1.90–1.88 Ga (Bergman et al., 2001; Sarlus et al., 2020; Andersson et al., 2021). During the initial phase of the Svecokarelian cycle, orogenic magmatism generated two comagmatic plutonic suites and associated supracrustal rocks with characteristic geochemical and petrological signatures (Bergman et al., 2001; Sarlus et al., 2019). Calc-alkaline granites, granodiorites, tonalites, and gabbro of the Haparanda Suite and equivalent volcanic rocks of the Porphyrite Group predominate in the eastern part of northern Norrbotten. In the western part, the alkali-calcic Perthite Monzonite Suite dominates with equivalent mafic to felsic volcanic and volcaniclastic rocks of the Kiirunavaara Group (porphyry group in Bergman et al., 2001). Related NE-SW–oriented crustal shortening (D₁ in Norrbotten, D₂ in the Skellefte district; Bergman and Weihed, 2020) caused basin inversion and a penetrative continuous tectonic foliation under plastic conditions (Andersson et al., 2020; Bergman and Weihed, 2020). The timing of the D₁ event was indicated to precede the Perthite Monzonite Suite based on radiometric dating of variably deformed plutonic rocks (Bergman et al., 2001), whereas it was interpreted from crosscutting relationships to be synchronous with the Perthite Monzonite Suite (Andersson, 2021). The following Svecokarelian magmatic processes during the second phase of the orogeny (ca. 1.82–1.78 Ga) generated I- to A-type granitic plutons of the Edefors Suite (as the early phase of the Transscandinavian igneous belt; Högdahl et al., 2004) and silica-rich S-type granites of the widespread Lina Suite (Öhlander and Skiöld, 1994; Bergman et al., 2006; Martinsson et al., 2016). The related D₂ deformation in Norrbotten is characterized by strong strain partitioning into reactivated shear zones under brittle-plastic conditions in response to E-W–oriented crustal shortening (Bergman et al., 2001; Andersson et al., 2021; Bauer et al., 2022). In the Gällivare area, a structural block situated between two crustal-scale shear zones (Nautanen deformation zone in the east: Bauer et al., 2022; and the Fjällåsen deformation zone in the west: Andersson et al., 2020) a gneissic S₁ cleavage was folded (F₂) during D₂, suggesting overprinting of the two deformation events (Bauer et al., 2018).

Dedicated studies on metamorphism from northern Norrbotten are insufficient to confidently assign peak events to the regional tectonic framework. Limited pressure-temperature (P-T) determinations utilizing mineral assemblages (chlorite-actinolite in ultramafic rocks, hornblende-plagioclase in mafic rocks) indicated a general transition from lower amphibolite facies in the west to upper amphibolite facies toward the east (Bergman et al., 2001). Based on mineral assemblages of mafic rocks, the Kiruna area (hornblende, actinolite, chlorite, plagioclase, albite, epidote) and the Stora Sjöfallet area (chlorite, actinolite, epidote) are isolated metamorphic domains with upper greenschist to lower amphibolite facies (Frietsch, 1979; Bergman et al., 2001; Skelton et al., 2018). Temporally, geochronological data of metamorphic ages indicate two regional metamorphic events (M₁ and M₂) that broadly overlap with the early and late Svecokarelian orogenic cycle (Bergman et al., 2001; Sarlus et al., 2018). The effect of contact metamorphism related to the magmatic suites is poorly constrained regionally, but in the Gällivare area, the timing of plutonism (1.81–1.78 Ga; Bergman and Weihed, 2020) and the amphibolite facies broadly correspond (Skelton et al., 2018; Sarlus et al., 2019). The regional hydrothermal alteration system is characterized by sodic-calcic-iron-chlorite metasomatism (Bergman et al., 2001; Martinsson et al., 2016; Logan et al., 2022). Alteration associated with IOA ore deposits is mainly ore-proximal sodic and distal potassic and is related to early orogenic events; however, a late potassic alteration can occur spatially with IOCG mineralization zones and D₂ shear zones (Martinsson et al., 2016; Bauer et al., 2022).

The Kiruna area is recognized as a prospective region for different mineralization styles, and economic deposits are variable in character, ranging from stratiform-stratabound and epigenetic-style sulfide deposits to IOCG-style, Kiruna-type (IOA), and skarn-rich iron ores (Martinsson, 1997; Martinsson et al., 2016; Logan et al., 2022). Hjalmar Lundbohm conducted the first geologic investigations in 1898 in the Kiruna district concerning the location, shape, and chemistry of the iron ores (Lundbohm, 1898). In the past century, the area has been extensively studied. Mapping campaigns (e.g., Ödman, 1957; Offerberg, 1967), stratigraphical/petrological studies (Frietsch, 1979; Martinsson, 1997), and geochronological (Cliff et al., 1990; Billström et al., 2019; Andersson et al., 2022) and structural investigations (e.g., Vollmer et al., 1984; Wright, 1988; Andersson et al., 2021) make the area relatively well understood.

Stratigraphy of the mining district: The depth to the Neoarchean basement rocks in the Kiruna area was approximated to be 2 to 3.5 km (Luth et al., 2018), increasing toward the south (Holmgren, 2013; Luth et al., 2018). A triple junction south of Kiruna led to north-northeastward propagation of the regional NNW-SSE–oriented rift and locally generated NE- to NNE-trending greenstone belts. This back-arc basin development in the Kiruna area provided space for the deposition of the Orosirian stratigraphic pile comprising poorly sorted conglomerates at the base, followed by bimodal volcanic-volcano-sedimentary rocks that transition into sedimentary basin-fill materials at the top (Andersson et al., 2021). In the Kiruna area, the greenstone belt is enclosed by younger Paleoproterozoic rocks and appears aligned to major shear zones (Bergman et al., 2001).

A thick pile of Rhyacian mafic to ultramafic volcanic rocks (2–4 km) constitutes the lowest part of the stratigraphic column (Fig. 2). At the base, the Kovo Group consists of calc-alkaline mafic to intermediate volcanic rocks and volcaniclastic sediments deposited on top of Archean rocks. The upper part of the Rhyacian, the Kiruna Greenstone Group is composed of six formations differentiated by petrographical and geochemical characteristics and contains mainly mafic to ultramafic metavolcanic rocks that formed during an intensive rifting phase (Martinsson, 1997). In general, the Kiruna Greenstone Group is characterized by initial clastic sedimentation and evaporite deposition that transitioned into extensive within plate basalt (WPB)-type volcanism. Outcrop observations indicate a northeast-southeast trend with a vertical to steep east dip (Fig. 3A). The stratiform-stratabound syngenetic Cu-(Zn) Viscaria deposit and the epigenetic Cu-Au Pahtohavare ores are situated within the Kiruna Greenstone Group (Martinsson, 1997; Logan et al., 2022).

Overlying the Kiruna greenstones, the Orosirian stratigraphic pile starts with poorly sorted, alluvial polymict conglomerates and graywackes of the Kurravaara conglomerate (Kumpulainen, 2000; Andersson et al., 2021), which were deposited at 1888 ± 3 Ma (minimum depositional age), marking the onset of Orosirian volcanism and basin development (Andersson et al., 2022). The Kurravaara conglomerate is followed upward in stratigraphy by the Kiirunavaara Group (Martinsson, 2004). The basal part of the Kiirunavaara Group constitutes trachyandesitic subvolcanic and extrusive rocks of the Hopukka Formation (also described as the “Syenite Porphyry” by Geijer, 1910) followed by rhyodacitic pyroclastic rocks of the Luossavaara Formation (also described by Geijer and Ödman, 1974, as the “quartz-bearing porphyry”). Outcrop observations in the northern part of the study area revealed undulating, bedding-parallel, NNW-SSE–oriented clasts within the Luossavaara Formation (Fig. 3B). At the contact between the Hopukka and Luossavaara Formations, the world-class Kiirunavaara IOA deposit forms a 5-km-long and approximately 100-m-thick, E-dipping (60°–70°) tabular magnetite-apatite body (Grip and Frietsch, 1973) that extends to approximately 2,400 m below ground level according to current knowledge (Luossavaara-Kiirunavaara AB, 2022). Thinner and shorter tabular orebodies are situated south (Konsuln, Sigrid, and Viktor) and north (Luossavaara) of the main mineralization at a similar stratigraphic position (Bergman et al., 2001; Luossavaara-Kiirunavaara AB, 2022). In the vicinity of the Luossavaara orebody, the Hopukka Formation contains ore-parallel magnetite dikes with oblique en echelon-style calcite veins (Fig. 3C). A syenite body (gabbro/monzonite in Bergman et al., 2001) situated parallel and stratigraphically below the orebody crops out within the Hopukka Formation, and in deeper parts, a granitic intrusion is emplaced (Westhues et al., 2016). Where observable in the underground mine, the contact between the granitic intrusion and its hanging wall is shallow E-dipping, and it remains open toward depth. The upper part of the Kiirunavaara Group consists of the Matojärvi Formation, which marks a change in the geologic development of the area and a transition to a volcanic-volcaniclastic-sedimentary sequence (Andersson et al., 2021). The volcanic portion of the Matojärvi Formation is bimodal and composed of basal rhyolitic tuffs followed by basaltic agglomerates, tuffs, and minor coherent volcanic rocks that transition into heterolithic breccia conglomerates, graywacke, and phyllites (Geijer, 1950; Martinsson, 2004; Andersson et al., 2021). Within the top of the Kiirunavaara Group, the Per Geijer hematite-magnetite-apatite–rich orebodies (Rektorn, Henry, Nukutus, Haukivaara, indicated in Fig. 1, and the newly discovered Per Geijer Deep situated underground) are situated at the contact between the Luossavaara and Matojärvi Formations (Andersson et al., 2021; Luossavaara-Kiirunavaara AB, 2022). The uppermost part of the Orosirian stratigraphy is marked by the Hauki quartzite, which is composed of cross-bedded arenitic rocks and horizons of breccia conglomerates (Martinsson, 2004; Andersson et al., 2021). Outcrop observations at the upper part of the Hauki quartzite indicate bedding-parallel foliation (Fig. 3D), and where bigger clasts are present from lower units, an undulating pattern with foliation wrapping around the clasts can be observed (Fig. 3E).

Characteristic alkali and calcic-iron alteration have spatiotemporal relationship to the IOA and IOCG mineralizations in the district (Martinsson et al., 2016; Andersson et al., 2021). The Kiirunavaara deposit and its host rocks have characteristic ore-proximal pervasive albitization that is crosscut by localized actinolite veins, and distal sodic-calcic (actinolite ± albite ± epidote ± titanite) to sodic-calcic-ferrous (actinolite and magnetite ± albite) alteration that transitions to potassic-ferroan (biotite ± magnetite) alteration (Martinsson et al., 2016; Paolillo and Giapis, 2021; Lupoli et al., 2022). The Per Geijer iron ores are mainly potassic altered (Westhues et al., 2016), with the host rocks containing occurrences of K-feldspar, sericite, chlorite, biotite, and carbonate (Martinsson et al., 2016). Ore-proximal alteration with a red-brick appearance from the Per Geijer area has been recorded as hematite inclusions within silica, affecting both footwall and hanging-wall rock types (Géhin and Miles, 2022). Regional metamorphism (M₁) is absent, and greenschist to lower amphibolite facies metamorphism was attributed to the late Svecokarelian events (M₂) in the Kiruna area (Bergman et al., 2001; Andersson et al., 2021).

Structural evolution of the Kiruna mining district: The area is characterized by uniform bedding attitudes over large distances with steeply E-dipping structural grain (Vollmer et al., 1984; Wright, 1988; Andersson et al., 2021). The stratigraphic pile contains a heterogeneously developed cleavage that is subparallel to the lithological contacts, and local shear zones together with highly tectonized units have a similar northeast to north-northeast strike (Andersson et al., 2021). The crustal-scale, NNE-trending, steep Kiruna-Naimakka deformation zone is adjacent to the Kiruna area (Bergman and Weihed, 2020).

Vollmer et al. (1984) suggested one major episode of crustal shortening (WNW-ESE) and argued that the Kiruna area is situated on the eastern limb of a major antiformal structure, with strain patterns influenced by diapirism of the neighboring granitic rocks. However, according to Wright (1988), intrusive contacts are undeformed, the cleavage is not consistent with the pluton contacts, and strain intensity cannot be correlated with the distance from the plutons. Regional structural mapping conducted by Wright (1988) distinguished four separate deformation events arguing for a fold-thrust model. Talbot and Koyi (1995) elaborated further on the fold-thrust model, arguing for the simultaneous conglomerate and quartzite deposition in a foreland basin. However, the Kiruna area, as well as the northern Norrbotten ore district as a whole, is characterized by steep deformation zones that are less likely to form owing to thrusting, and evidence for simultaneous rotation of crustal-scale faults is lacking both in the Kiruna mining district and regionally (Bergman et al., 2001). The latest established structural framework model (Fig. 4) by Andersson et al. (2021) details the structural evolution of the Kiruna mining district in accordance with four deformational events that are in line with regional tectonic events as well as the petrological character of northern Norrbotten. Back-arc extension and basin development synchronous to emplacement of IOA mineralization (D₀) was followed by basin inversion with crustal shortening in the east-west direction (local D₁, regional D₂) and refolding in response to north-south crustal shortening (D₃) and later (D₄ related) hydraulic fracturing (Andersson et al., 2021). The Kiirunavaara Group is lacking the regional S₁ fabric, or it was overprinted within the Kiruna mining district area (Andersson, 2021). However, based on kinematic evidence, N-S–trending shear zones were active west of Kiruna during D₁ (Andersson et al., 2020), and folded S₁ fabrics within the Kiruna Greenstone Group have been identified 5 km south of the Kiirunavaara IOA deposit (Logan et al., 2023), indicating differential fabric development within the stratigraphic succession of the mining district.

Basin inversion during D₂ crustal shortening (E-W to NW-SE directed) reactivated preexisting normal faults and generated moderate to steep NNE-SSW–trending reverse oblique to dip-slip brittle-ductile high-strain zones, developed at lithological contacts and in units susceptible to elastic-plastic behavior (Kurravaara conglomerate, Matojärvi Formation, and the Hauki quartzite). The wide range of rock types caused strong competency differences, leading to strain partitioning during the D₂ crustal shortening. The IOA orebodies situated at lithological contacts acted as rigid bodies, and they were interpreted to be lineation-parallel boudins (Vollmer et al., 1984; Wright, 1988), similar to mesoscale boudinage and pinch-and-swell structures of quartz veins from the Kiruna area (Andersson et al., 2021). Pinch-and-swell structures of quartz veins observed in the southeast of the area within intrusive Perthite Monzonite Suite granites (Fig. 3F) indicate that D₂-related compression affected the intrusive rocks as well, but the extent is unknown. Volcanic and volcaniclastic rocks are intensely foliated within and along the margins of the rigid ores, suggesting a high-strain gradient (Wright, 1988). Within competent volcanic units and orebodies, the crustal shortening was mainly accommodated by brittle fracture planes. These spatiotemporal associations of the lithological units during structural evolution, expressed by the rheological heterogeneity and the interplay between ductile and brittle structures during different deformation events, played an important role in controlling the crustal architecture (Andersson et al., 2021).

To characterize the subsurface, geologists recognized early the need to view the world as multidimensional (Turner et al., 2020), and fast technological advancements continuously favor this transition. The geologic features of a mining district can be represented using a three-dimensional geologic model that integrates a multitude of surface and subsurface data such as geologic, geophysical, geochemical, and geotechnical measurements and observations. Geologic rules and knowledge can be integrated into geologic models through the consideration of conceptual geologic models (Unger et al., 2018). The importance of conceptual models and geophysical data is proportional to the increasing distance away from the densely drilled areas, especially in regions with poor outcrop exposure.

Currently, geologic modeling is based mainly on explicit and implicit modeling methods (Cowan et al., 2002; Caumon et al., 2009; Stoch et al., 2018). As with computer-aided design (CAD) applications, traditional hand-drawn cross sections and surface maps can be transferred into the digital world by a manual (explicit) drawing approach. The interpreted geologic features from parallel sections can be connected to generate 3-D triangulated surfaces (Mallet, 2004; Caumon et al., 2009) facilitating the transition from 2-D to 3-D. These 3-D discrete elements can enhance insights into the overall structural architecture of an area. However, they are inherently subjective, nonunique, and time-consuming owing to the manual input and the high uncertainty between parallel cross sections or data points, leading to simplifications and low model resolutions (Schetselaar, 2013; Stoch et al., 2018). The nonreproducibility of the explicit models is related to the hermeneutic components of geology, which makes it highly interpretative, as observations cannot be separated from geologic knowledge (Wellmann and Caumon, 2018). Implicit modeling techniques aim to mathematically describe the geologic surfaces in 3-D space by computing a scalar field, where a geologic surface location will have the same scalar value. Multiple surfaces can be represented as the isovalues of a single or several scalar fields (Lajaunie et al., 1997; Wellmann and Caumon, 2018). Implicit surfaces can be created using meshless methods such as radial basis function (RBF) interpolation (Carr et al., 2001; Cowan et al., 2003; Hillier et al., 2014) or dual kriging (Lajaunie et al., 1997; de la Varga and Wellmann, 2016) and mesh-based methods (Moyen et al., 2004; Caumon and Collon-Drouaillet, 2014). These techniques permit the building of geologic models in areas with a high degree of complexity (Wellmann et al., 2017), making them reproducible and less biased (Grose et al., 2019). Geologic observations and measurements are incorporated as spatial data (such as points, polylines, or structural data) and are fitted to the implicit function (Vollgger et al., 2015). Planar structural data can be incorporated as the gradient of the scalar field at specific locations (Hillier et al., 2013), or structural frames can be employed to characterize the geometry of major structural features and used in time-aware geologic modeling (Grose et al., 2021a, b). Several commercial software (e.g., Leapfrog by Seequent, Gocad-SKUA by AspenTech, and 3-D GeoModeller by the French Geological Survey and Intrepid Geophysics) and open-source packages (e.g., Loop-Structural presented in Grose et al., 2021a, b; GemPy presented in de la Varga et al., 2019) exist, and they employ a diverse range of algorithms and workflows for interpolation.

In mining districts, data-rich environments such as the vicinity of the mined orebodies can provide substantial control points, whereas other areas are underconstrained owing to targeted exploration and exploitation of mineral resources. Creating a geologic model of a mining district with heterogeneous data distribution makes it necessary to include geologic rules and experience by using rule-based or conditional implicit modeling techniques (Unger et al., 2018). In this way, the 3-D model is created by a surface-based modeling approach, where the interpolated unique contacts are lithological boundaries (Stoch et al., 2018) describing the varying scalar physical quantities in space (Carr et al., 2001; Cowan et al., 2002, 2003). Implicit techniques can generate models that predict the geologic realm to a higher degree in areas with sparse data sets by the addition of synthetic data (i.e., additional constraints that are not backed by direct observation) based on the geologic knowledge of the modeler (Jessell et al., 2014).

The geologic model of the Kiruna mining district was done using the Leapfrog Geo 2022.1 software package, which is deemed suitable for modeling areas with high spatial data heterogeneity. The software package uses radial basis functions (RBFs) to build 3-D geologic models by spatially interpolating the data set (Cowan et al., 2003; Stewart et al., 2014). The Applied Research Associates, New Zealand (ARANZ), in collaboration with the University of Canterbury, developed the FastRBF, which made possible the interpolation of larger data sets characteristic to geologic studies. ARANZ (now Seequent) began utilizing the basis function (the spheroidal function), and it was integrated into the Leapfrog software (Beatson et al., 1999; Carr et al., 2001; Stewart et al., 2014). This function has linear behaviors near the origin, whereas it has an asymptotic behavior away from the data, contrasting other basis functions that asymptotically approach a sill value and have a smooth form near the origin (Stewart et al., 2014).

The data set used for modeling includes surface topography, geologic map, aeromagnetic anomaly map, underground mine wall mapping, drill core logs, and outcrop structural measurements. The volume of the geologic model (Fig. 5) stretches over 11.7 km (northing) by 9.5 km (easting) by 3 km (depth), encompassing all the geologic units from the Kiruna area. The surface topography was triangulated based on Lantmäteriet elevation data (2+-m grid precision) and shaped to the extent of the model area using mesh clipping operations. The underground mine wall mapping data and drill core logs containing lithological, structural, and geotechnical data were supplied by Luossavaara-Kiirunavaara AB (LKAB). The area of focus can be subdivided based on the availability of data into two densely drilled parts (Fig. 5), consisting of ~6,000 drill holes (total drilling length of ~1,100 km): the Kiirunavaara area, which is currently the only active iron ore mine in the district, and the Per Geijer deposits, situated a few kilometers farther north. From the ongoing underground mining operation at the Kiirunavaara mine, underground mine wall mapping data contained ~92,000 measurements of ore contact, joint, and joint fill (Fig. 5). Outside the two data-rich areas, current geologic knowledge is based mainly on outcrop observations and a few drill cores available from the Swedish Geological Survey (SGU). The subsurface measurements and observations were imported, upscaled, and used as primary data (App. Table A1), whereas synthetic data (App. Table A2) was created to guide the lithological boundaries and faults, according to field observations and established conceptual models (Offerberg, 1967; Vollmer et al., 1984; Wright, 1988; Martinsson and Perdahl, 1993; Talbot and Koyi, 1995; Bergman et al., 2001; Andersson et al., 2021). The synthetic data served as additional interpretative points in locations where subsurface data was not available, and it was added as structural discs specifying the location and orientation of surfaces.

The geologic map presented in Figure 1 was used for modeling, and it was modified based on the work of Offerberg (1967), Bergman et al. (2001), and Andersson et al. (2021) and supplemented by our own field observations and measurements. In locations where the lithological contacts on existing geologic maps were not consistent, the boundaries were inferred from the aeromagnetic anomaly map (App. Fig. A1) and used together with field observations of the lithostratigraphic units (Fig. 3) to guide the addition of synthetic data.

The reduction of spatial constraints in areas situated outside of the densely drilled areas and outcrop observations increases the uncertainty of the model. The proximity to drill holes and outcrop observations is presented as distance function contour lines in two sections A-Aʹ and B-Bʹ (Fig. 5) and highlights the areas where the models rely on geophysical and conceptual models. The distance function can be considered a geodiversity metric (Lindsay et al., 2013) that evaluates the contribution of data to model variability, and the spatial data density maps (Fig. 5) pinpoint areas where data collection must be prioritized in future geologic mapping campaigns. Lindsay et al. (2013) indicated that multiple model realizations permit uncertainty simulation and lead to a better understanding of the geometric variability and target geology. The synthetic data set in the current study aims to provide the best possible set of orientation measurements according to the current knowledge. As a range of other scenarios might exist that better represent the subsurface, the synthetic data set can be substituted and a suite of other geologic models can be generated, without modifying the well-constrained areas.

De Kemp et al. (2016) highlight the necessity of a workflow in geologic modeling that can help to focus on the defined tasks and assure data consistency and purpose-driven geologic modeling. The presented workflow highlights the general modeling procedure and serves as a guide in achieving a robust geologic 3-D model that is consistent with the input data and the geologic framework. The workflow that was followed during the modeling process is presented in Figure 6. Instead of a linear approach, the workflow is based on a lateral approach where geologic knowledge has a higher weight, and the model can be restructured and changed as the underlying conceptual model progresses. The final model, which is situated at the core of the process, will be affected by the change of any peripheral component of the workflow. This modeling methodology proved to be appropriate for active mining areas where new exploration campaigns can bring new data and the model can be updated periodically and for geologically complex areas where established conceptual models are frequently challenged by new models. Variably constraining the model elements, instead of modeling only the parts from densely drilled areas, provides a valuable addition by enhancing the overall structural architecture and optimizing the integral perception of the 3-D geologic models.

1. The first two steps included the import, examination, and standardization of the collected data. The geologic map, drill hole, and underground and surface measurements were converted to a common coordinate system—the Swedish national grid SWEREF99 with local projection to the local meridian of 20° 15″. In instances where missing values were detected for a location or drill hole path data, the affected holes or mine wall mapping measurements were excluded from the study. The proportion of data discarded at this stage was negligible, constituting less than 1% and primarily attributed to historical drill hole data sets.

2. The third step includes data visualization in 3-D, which was used to identify and mitigate errors due to incorrect location, duplicates, or incorrect georeferencing. This aimed at confirming the alignment of the geologic map, topography, and satellite orthophotograph and further confirmed the accurate placement of surface measurements and drill hole collar locations against established control points such as outcrop locations, lakes, and transportation networks. The relative location of mine wall mapping measurement locations (Fig. 5) was inspected, and 137 locations that were outside of the mine area were removed. By geometrical pattern analysis of the drill core logs using the 3-D viewer (similar to Cowan, 2020), the distribution and trend of the data were visually analyzed. Data from inconsistent drill core logs (e.g., where historical lithological logs are not consistent with adjacent and new logs) or data with limited coverage due to modifications in logging procedures (e.g., measurements and observations encompassing only a partial area resulting in gaps in certain sections) were excluded from the analysis and not considered for modeling. The subsequent analysis focused on systematically recorded lithological, structural, and geotechnical logs spanning the entire study area.

3. During step four, a representative fault network was interpreted from geologic maps and aeromagnetic anomaly maps and assessed in the context of conceptual structural framework models described in the literature (cf. Parák, 1969; Vollmer et al., 1984; Wright, 1988; Andersson et al., 2021). The geometry of the geologic structures was traced by stratigraphic evidence, structural markers, and structural measurements from drill holes, whereas on the surface it was traced from structural measurements and observations in outcrops. The major deformation zones are developed within NNE-striking lithological contacts, and these shear zones overlap with lithological boundaries, making it possible to use stratigraphic contacts from drill core logs as geometry indicators together with markers indicative of shearing (e.g., mylonite). Brittle NE- and NW-striking faults were delineated based on drill core makers with an indication of faulting, core loss, rock quality description (RQD) values, and (in certain cases) breccias. In the case of brittle faults, several markers indicative of faulting were observed with a consistent orientation that can be traced throughout several drill holes. The faults were traced on the topographic surface mesh, whereas the three-dimensional attitude of the fault cores was assigned as structural trends based on structural analysis of the areas with frequent fault markers. An example of a zone with 76 drill core fault markers from a cluster of WNW-ESE–trending faults situated between the Rektorn and Henry open pits is presented in Figure 5. The stratigraphic column (Fig. 2) was modified after Martinsson et al. (1999) and Andersson et al. (2021), and it was used to define crosscutting and chronological relationships of the lithostratigraphic units. To establish a contact between two units, unique and unambiguous contact surfaces are needed; however, lithological logs usually contain the rock types. In the Kiruna mining district, the contacts between units are not always obvious, and the resultant model cannot represent the complexity of very detailed logging. This required lithological codes from logging to be grouped into the five main units that are intercepted by drill holes. Given the well-preserved stratigraphic succession in the Kiruna area and the established lithological column, the logged lithological codes align effectively with the stratigraphic column (Fig. 2). Grouped rock types show a typical “barcoding” of alternating layers of the different volcanic and sedimentary units, and intervals less than 1 m were not considered. Units that are not intercepted by drill cores were constrained using the geologic map and surface outcrop measurements and further guided by the addition of synthetic data.

4. The geologic model was built during step five using the refined data and established geologic knowledge. Nine E- to SE-dipping supracrustal units were modeled in three dimensions, and they define the internal lithostratigraphic volumes. Structural repetitions are modeled as individual lithostratigraphic volumes, and together with intrusive units a total of 26 3-D elements were created. The NNE-striking structures developed frequently at lithological contacts and certain formations represent highly strained zones. These shear zones were not modeled owing to their extensive nature, and current software packages are not capable of factoring in the competency difference of rock types or areas affected by polyphase deformation events. In total, nine NW-SE– and NE-SW–trending, steeply dipping faults were constructed and incorporated into the fault network, considering their crosscutting and terminating relationships. The slip direction of the faults is uncertain, and they were modeled as static explicit faults. The geologic formations were modeled based on the lithological and structural logs and contact relationships (deposition, erosion, intrusion), and younging directions of the surfaces were defined by ordering them in chronological order according to the stratigraphic column (Fig. 2). Trend directions (App. Table A3) were applied based on the E-dipping angles of the structural grain of the area that generated an ellipsoidal anisotropy. Using the surfaces from the lithological boundaries and additional synthetic data, a lithological boundary model of the formations was created using the hybrid explicit-implicit modeling method.

Ongoing mining operations and recent exploration efforts within the Kiruna mining district led to the collection of significant data sets from the surface and at various depths, offering heterogeneously distributed spatial data from the mining district. Deposit-scale structural analysis and interpretation of geotechnical log parameters aim at complementing previous surface-based studies and highlighting the 3-D relationships between the elements of the crustal architecture. In the following section, the structural characteristics of the local stratigraphy will be presented from subsurface data, and the structural analysis is divided into two parts: the first area includes the stratigraphically lower Kiirunavaara deposit and its host rocks (Hopukka-Luossavaara Formations), whereas the second area is situated higher up in the stratigraphy and encompasses the Per Geijer deposits and associated formations (Luossavaara-Matojärvi Formations). These parts are further subdivided into domains separated by major brittle faults (Fig. 5). Structural measurements were analyzed using the Leapfrog Geo stereonet analysis tool (equal-area, lower-hemisphere).

Kiirunavaara area: Underground mine wall mapping from the Kiirunavaara mine constitutes the most abundant data set from the area, and measurements of fractures can be separated into three main sets with different orientations (Fig. 7A). Generally, the three predominant fracture sets include a moderately steep SSW-dipping (108°/60°) structure and two approximately N-S–trending, east (351°/70°)- and west (199°/75°)-dipping structures (Fig. 7A). In outcrops west of the Kiirunavaara mine, fracture sets with the same three orientations were observed in basaltic rocks of the Kiruna Greenstone Group (App. Fig. A2A). In addition to the main group of fractures, several subsidiary sets can be observed, trending northwest-southeast and dipping both northeast and southwest. The SSW-dipping fractures (set 1 in Fig. 7A) are numerically dominant throughout the mine area and contain mineral fills predominantly of anhydrite, chlorite, calcite, feldspar, pyrite, clay minerals, and less abundant hematite and biotite (Fig. 7B; App. Fig A2B). The approximately N-S–trending, W-dipping fractures (set 3 in Fig. 7A) are the second most abundant, and together with the E-dipping set (set 2 in Fig. 7A), they contain mineral fills of mainly hematite, biotite, and chlorite (Fig. 7B; App. Fig. A2C, D). The spatial relationship of the fracture sets has been examined based on four fault domains (Fig. 7C) separated by WNW-ESE– and WSW-ENE–trending brittle faults that in some cases crosscut the Kiirunavaara orebody. The main fracture sets have distinct orientations in individual fault domains throughout the Kiirunavaara area that can be followed using characteristic mineral fills. Set one has a SSE-dipping character in the southern part, which gradually shifts toward the north toward a south-southwest dip. Set two changes from a northeast dip in the south toward a steeply E-dipping character in the north. Fracture set three appears mainly in the northern part, dipping west, and its frequency gradually decreases toward the south. Within the Konsuln satellite orebody, fracture set three appears as a weak ESE-dipping structure.

Structural measurements from oriented drill cores are available mainly from the northern end of the iron ore deposit and the hanging-wall rocks, reaching north toward the Luossavaara deposit. Measurements of foliation, banding, and ore contacts follow a similar trend, with N-S–trending and moderately steep (45°–65°) E-dipping structures throughout the deposit (Fig. 8A). Two fracture sets are present, set one with a moderately steep east-southeast dip, and a second steep E-dipping set with decreasing frequency toward the north (Fig. 8A). The mineral fill of the first fracture set (App. Fig. A2C) is dominated by anhydrite, hematite, pyrite, chlorite, quartz, and talc with actinolite and calcite as subordinate minerals, whereas the second fracture set (App. Fig. A2B, D) contains mainly chlorite, actinolite, calcite, and hematite (Fig. 8B). Measurements of lineation on fracture planes (Fig. 8C) are measured as gamma values from the oriented drill cores. Lineations on the first fracture set have a shallow plunge toward the east-southeast and southwest, whereas the second fracture set is characterized by steep to moderately steep plunges toward the south or less frequently toward the north.

Rock quality description and core orientation quality based on continuous intervals that were fitted together from core runs are presented in Figure 9A and B for the Kiirunavaara area. Based on three E-W–oriented cross sections through the orebody, drill core markers with low orientation quality (discontinuous orientation between drilling runs) are situated close to the ore contact. Markers with high orientation quality (continuous orientation between drilling runs) are situated mainly within the orebody or its host rocks (Fig. 9A). On the same cross sections, RQD values show similar results, with very poor (0–25%) values situated mainly at the ore contact, whereas a gradual increase can be observed toward the hanging wall and footwall (Fig. 9B).

Per Geijer area: Structural measurements from the oriented drill cores (Fig. 10A) have been separated into three different domains with a similar approach as that for for the Kiirunavaara deposit. The measurements are taken mainly from the proximity of the Hauki, Rektorn, Henry, and Nukutus open pits and therefore are characteristic only for the Hauki quartzite, Matojärvi Formation, and to a lesser extent the Luossavaara Formation where the drilling stops. Planar structural measurements of foliation, lamination, veins, banding, and lithological contacts follow the same north-south trend and dip moderately steeply toward the east. Close to the Henry and Nukutus open pits, measurements indicate a steepening and a gentle shift of the strike direction toward the northeast, which is in alignment with the general strike of the lithological boundaries on the geologic map (Fig. 1). Structural measurements were recorded from markers indicating both ductile and brittle deformation. Drill core markers indicating brittle faulting constitute two main sets (Fig. 10A). The first set defines a continuous trend (north of Rektorn) that can be traced across several drill holes, having a northwest-southeast trend with a steep northeast dip. The second set has a consistent north-south trend and a moderately steep dip toward the east throughout the area. Drill core markers indicate that mylonitic fabrics are trending north-south and dip steeply toward the east.

Rock quality descriptions from the Per Geijer area (Fig. 10B) show a heterogeneous spatial distribution compared to the Kiirunavaara area, and drill core markers with low values dominate throughout the Hauki quartzite and Matojärvi Formation. Characteristic subareas indicative of strong or weak rock mass that cannot be distinguished as low RQD values are frequent and occur homogeneously, associated with a limited number of high RQD values. The elevated degree of fracturing and shearing of the rock mass is characteristic of the whole area.

The following section describes the structural framework and lithological boundary 3-D models of the Kiruna mining district, focusing on the stratigraphic and tectonic relationship between the lithostratigraphic units and structural elements. The elements of the structural framework (Fig. 11A) include the dominant faults and shear zones that have a major impact on the geometry of the lithological units and orebodies. Markers with indications of shearing (mylonite, core loss, and to a lesser extent low RQD values) present continuity in N-S–trending directions throughout the study area. These indicators have been recorded frequently at the lithological contacts, indicated by low orientation quality and RQD values in the Kiirunavaara area (Fig. 9), or distributed homogeneously within less competent units such as Matojärvi Formation, indicated by dispersed low RQD measurements (Fig. 10B). The N-S–trending shear zones have not been modeled as fault planes owing to their wide extent and close relationship with the lithological contacts. However, the boundary surfaces of the lithological contacts are included in the structural framework model (Fig. 11A). These N-S–trending shear zones have been described as synvolcanic normal faults developed during basin formation that experienced reverse reactivation during crustal shortening (Andersson et al., 2001). The brittle fault network includes five NW-SE– and four SW-NE–trending faults that crosscut each other at the extremity of the orebodies, creating complex conjugate structures. At the southern end of the Kiirunavaara orebody (Fig. 11A) the faults lack structural constraints from the subsurface, but abrupt lithological changes can be observed in lithological logs, and they constitute clear geophysical lineaments displacing the syenite intrusion and the Konsuln orebody relative to the main intrusion and the Kiirunavaara orebody (Fig. 11A). At the northern end of the Kiirunavaara orebody, the presence of drill core markers indicative of faulting increases, and abrupt changes in the lithological logs occur, highlighting the presence of faults in the area. The SW-NE–trending faults of the conjugate set are dipping toward the southeast and are located close to the southern tip of the Luossavaara orebody. The NW-SE–trending fault of the same conjugate set is dipping steeply toward the northeast. The shape of the Kiirunavaara orebody does not follow the fault orientation, but above the northern end, a continuous band of hematite can be traced through several drill holes near massive magnetite, potentially indicating the presence of a fault or a shear zone. In the Per Geijer area, between the Rektorn and Henry open pits from the Per Geijer area, the WNW-ESE–trending fault (Fig. 10A) is well constrained from subsurface markers and has orientations similar to those of NW-SE–trending faults from the conjugate sets in the south of the study area; however, a southwest-northeast component has not been identified.

Nine supracrustal units were modeled in three dimensions and the volumetric lithological boundary model is presented in Figure 11 as an exploded-view diagram showing the relationship and assembly of various parts. The geologic model is made up of eight intrusive, 14 supracrustal lithological volumes, and four larger IOA horizons. The general attitude of the lithological units is north-south trending and steep to moderately steep east dipping, but local deviations occur, which are evident from the geologic map (Fig. 1) and aeromagnetic map (App. Fig. A1). From the modeled volumes (Fig. 12) it can be observed that the geometry of the basin changes upward in stratigraphy. The Kiruna Greenstone Group defines a NNE-trending, steeply ESE-dipping succession, whereas after the Hopukka Formation, which forms the largest unit by volume, the attitude of the lithological units changes to north trending, moderately steep east dipping. Supraposed Rhyacian basaltic units are north trending and steeply east dipping. The Kiirunavaara, Luossavaara, and Konsuln orebodies situated at the contact between the Hopukka and Luossavaara Formations constitute individual orebodies, and they are separated by sets of conjugate faults. In the proximity of the conjugate faults at the northern end of the Kiirunavaara orebody at approximately –1,700-m depth, the lithological contacts show a small-scale structural repetition of the Luossavaara Formation, where the rhyodacitic rocks are thrusted under the orebody and locally become the footwall to the ore. Two perpendicular cross sections of the geologic model illustrating the structural repetition are presented in Figure 13. The Matojärvi Formation is the smallest supracrustal unit by volume, is primarily found in the Kiruna area, and is absent regionally (Bergman et al., 2001). Despite its limited size and regional presence, the Matojärvi Formation represents the most heterogeneous unit. The Per Geijer ores are situated within or at the lithological contact of the Luossavaara and the Matojärvi Formation, forming several IOA orebodies of different sizes.

Over the past century, the geologic setting of the area, including stratigraphy (Geijer, 1919; Parák, 1975; Frietsch, 1979; Martinsson, 1997) and the structural framework (Vollmer et al., 1984; Wright, 1988; Talbot and Koyi, 1995; Andersson et al., 2021, 2022) have been intensely studied and resulted in several contrasting interpretations. Previous studies within the mining district were based predominantly on surface mapping (e.g., Vollmer et al., 1984; Wright, 1988; Andersson et al., 2021) or to a lesser extent on mine wall mapping (Berglund and Andersson, 2013). The 3-D models presented in the current study are based predominantly on data obtained from drill cores and provide an important vertical range for the geology of the Kiruna mining district, making it possible to test previous and future conceptual geologic models of the area. This is of particular importance since the mining district is accepted as the key locality for IOA deposits, and the 3-D models can be used as a base when multidisciplinary mineral system modeling is conducted and when different genetic models are evaluated.

The geologic models of the Kiruna mining district presented in this study (Figs. 11A, 12) offer a 3-D insight into the crustal architecture and spatial association between the IOA orebodies, lithostratigraphic units, and brittle to ductile high-strain zones (Fig. 11B). Although 2-D geologic and geophysical maps provide foundational data for 3-D models, they are limited in offering a comprehensive insight into the subsurface, and a 3-D geologic model that synthesizes various data sets can enhance the understanding of the crustal configuration. The heterogeneity in data resolution within the mining district generates areas with sparse data distribution, which makes it difficult to apply exclusively implicit or exclusively explicit geologic modeling methods. A rule-based (Stoch et al., 2018) hybrid modeling technique was applied to model the Kiruna mining district. This modeling approach proved to be efficient in the present study in which the input data was continuously collected throughout the modeling process, and owing to ongoing exploration activities, additional data imported into the model allows for continuous updates and refined interpolants. After standardization and inspection of the data sets from densely sampled areas, the resultant geometries of the surfaces created by implicit geologic modeling represent the subsurface with greater detail. The interpolation of lithological logs or structural measurements in areas with distant data spacing by implicit modeling techniques generates surfaces that are not representative and often geologically incorrect. The incorporation of additional synthetic data by modeling geologists can help in generating geologically acceptable surfaces that are consistent with the established conceptual models. During the creation of geologic models of a mining district, the collected data and the synthetic data can be treated as separate data sets but used to model the same continuous surface. The current study presents one model realization for the Kiruna mining district, and a potential limitation of the hybrid implicit-explicit method is related to the time required to produce several model realizations. However, the method facilitates the possibility to test multiple geologic scenarios by replacing the synthetic data set in the future as the geologic knowledge expands. Using this approach, regional geologic models can be created and local high-resolution models of mining areas can be connected. Local models embody the descriptive characteristics of the larger geologic environment, and the integration of geologic models with different resolutions can improve the understanding of the regional geologic context.

Lithological contacts in the Kiruna mining district are not always sharp; therefore, characterizing the transition between lithofacies classes with gradational transition (i.e., alternating layers) is not representative in many cases and requires simplification of the logging results or different modeling techniques. Investigating lithostratigraphic units with high internal heterogeneity and transitional contacts could lead to better results by modeling the mineral abundance or creating lithological boundaries with buffer zones that are characteristic of the alternating layers. Integrated investigations including geochemistry, mineral compositional data, textural information, and petrophysical properties with existing geologic models can provide a framework to integrate geologic processes within geophysical surveys. The current geologic models of the Kiruna area serve as a starting model for future integration with multidisciplinary data sets and regional models, which can help in improving the knowledge of the mineral system and aid the search for new mineralization zones.

The eastward dip of the NE-SW– to N-S–trending lithostratigraphic units and shear zones within the Kiruna area are generally decreasing upward in stratigraphy (eastward on the geologic map) and are bound by a major stratigraphic repetition to the east (Fig. 12). A fold-thrust belt model involving WNW-directed crustal shortening was proposed by Wright (1988) and further developed by Talbot and Koyi (1995) but was rejected by Bergman et al. (2001) because of the steep nature of the shear zones. The structural framework and lithological boundary models (Fig. 12) indicate that the well-constrained contact between the Hopukka and Luossavaara Formations is dipping 55° to 65° toward the east. Lower angles can be observed only in the case of the Hauki quartzite, which varies between 50° and 70°, dipping toward the east. Steep angles inferred from the geologic model support previous studies suggesting reactivated normal faults (Andersson et al., 2021) associated with extensional setting during back-arc basin development (Sarlus et al., 2018; Andersson et al., 2021).

It was suggested by Vollmer et al. (1984) that the Kiirunavaara orebody and proximal satellite orebodies once formed a continuous horizon that was boudinaged during ductile deformation, and Andersson et al. (2021) indicate that lineation-parallel boudins served as rigid bodies during deformation. Wright (1988) described boudinaged magnetite veins in the proximity of the Luossavaara deposit, indicating that individual segments are separated and that the absence of necking is a result of high ductility contrast. In the central Kiruna area, early Svecokarelian-related D₁ fabrics are not present or masked (cf. Andersson et al., 2021; Logan et al., 2022), and the area has been suggested to represent shallow crustal levels where plastic deformation was not recorded, whereas late Svecokarelian-related D₂ occurred under brittle-ductile conditions (Andersson et al., 2021). The presence of ~E-W–oriented faults has been indicated by Parák (1975), suggesting that they might control the shape and position of some ore lenses; however, descriptions of these faults were not presented. The structural framework model (Fig. 11A) illustrates the approximately E-W–trending fault surfaces constrained by drill core markers, structural measurements from oriented drill cores, and aeromagnetic geophysical maps, highlighting the relationship between the ore lenses and the conjugate faults that separate these orebodies and create individual fault blocks. Furthermore, the topological relationship between the principal N-S–striking high-strain zones and the conjugate fault networks is evident in the structural framework model (Fig. 11A), and their proximity to the IOA ore lenses suggests that they might have influenced the position of the present-day IOA ore lenses. Steep to moderately steep SW-, SE- and S-plunging lineation measurements (Fig. 8C) indicate that movement occurred across faults with varying orientations.

Based on surface observations, Andersson et al. (2021) suggested the presence of Riedel shear fractures that developed in competent volcanic rocks in the Kiruna area as a result of E-W–oriented crustal shortening under brittle-ductile conditions with reverse, east-side-up sense of shear. The three dominant fracture sets observed within the Kiirunavaara orebody (NNW-SSE, NNE-SSW, and WNW-ESE) mimic a Riedel system (Fig. 7), implying that the Kiirunavaara orebody was fractured during the D₂ basin inversion event. However, as geochronological U-Pb titanite data indicate some brittle ore-parallel structures to be synchronous with the ore emplacement (Andersson et al., 2022), it cannot be excluded that parts of the fault network were formed much earlier and experienced tectonic reactivations during the regional D₂ event. South-plunging lineation measurements are more abundant, indicating that displacement on ore-parallel fractures was stronger, whereas hydrothermal mineral fill of anhydrite, pyrite, chlorite, and clay minerals occurs more frequently within WNW-ESE–trending fractures, suggesting that these fractures channeled fluids to a larger extent. This is further supported by the decrease in drill core orientation quality and RQD at the Kiirunavaara ore contact, suggesting the presence of distinct ore-parallel structures. Fractures both parallel and at high angles to the Kiirunavaara ore have been reported from a local-scale study of the mine, and the largest displacement was observed on ore-parallel structures (Berglund and Andersson, 2013). The slight variation in the angular relationship between fracture sets situated in different fault blocks (Fig. 7C) indicates different strain accommodation across the subdomains of the orebodies. These subdomains are separated by ~E-W–trending conjugate brittle faults that could accommodate the movement between individual blocks. Both oblique slip and dip slip were observed along similar latitudes in the proximity of the Luossavaara IOA deposit, and the slight difference between the reverse kinematics recorded by the structures was attributed to either stain partitioning or cyclic reactivation of different parts of the system (Andersson et al., 2021). The local structural repetition (Fig. 13), where the hanging-wall rocks are thrusted under the orebody, is situated at the intersection of a conjugate fault set, south of the Luossavaara ore (Fig. 11A). The thrust might be the result of an interplay of ore-parallel, oblique to dip-slip shear zones and additional vertical movement along NW-SE– and SW-NE–trending fault systems. Nevertheless, the small extent of the structural repetition suggests that limited movement was accommodated along the fault planes. The position of the Kiirunavaara ore and its satellite orebodies (Luossavaara and Konsuln) in map view combined with the lineation orientations suggest that synchronous displacement in the conjugate fault networks is kinematically compatible with the separation of the boudins and juxtapositions of individual blocks from slightly different crustal depths (Fig. 11B).

The implications of this study highlight the role of 3-D geologic modeling in elucidating the district structural and lithological framework. Integration of subsurface and surface data sets by a hybrid implicit-explicit geologic modeling method offers insights into the larger geologic environment and local mineralization zones. The close spatial relationship between N-S–trending shear zones, ~E-W–trending conjugate fault networks, and the IOA orebodies inferred from the geologic models and structural analysis provides insights into the structural controls and spatial arrangements of the mineralization. Nevertheless, to create a full 2-D or 3-D balanced restoration, shear strain and strain accumulation in different rock types must be further investigated, as strain partitioning has been a significant factor influencing the deformation processes.

The results of this study include three-dimensional lithological and structural framework models that integrate extensive mapping and drill hole data sets from the surface and subsurface. The models display the internal structure of the supracrustal rocks and reveal the spatial relationship between N-S–trending shear zones and NW-SE– to SW-NE–striking conjugate fault networks that control the present-day location of IOA ore lenses. A rule-based hybrid implicit-explicit modeling technique has been used to fully constrain the models with the existing data and the conceptual framework. Using this approach, it is possible to define geologic conditions like surface chronology, contacting relationships, and structural trends. In areas where drill holes or surface measurements were not available, geologic knowledge was implemented using synthetic data. Using this new approach, the artificial data set can be isolated from the rest of the data and is used only to implement the geologic knowledge where a lack of measurement and observations exists. As exploration advances and geologic knowledge is improved, synthetic data can be modified or replaced by real measurements and observations. In areas where direct geologic observations are limited, several conceptual and structural framework models can exist for the same area, and they can be separately implemented by replacing only the synthetic data set, producing several geologically realistic models for exploration and mining purposes. Parts of the model that are created by implicit methods using the densely sampled, measured data sets will remain consistent and unaffected. The presented geologic models provide valuable insights into the relationship between lithological units and structures, highlighting the proximity of conjugate fault networks with the IOA ore lenses. Deposit-scale structural analysis based on subsurface data indicates that competency contrast between heterogeneous lithological units played an important role during basin inversion. In accordance with previous surface-based studies, drill core markers suggest that incompetent domains such as the Matojärvi Formation accommodated noncoaxial strain to a higher extent compared to competent domains such as the iron ore lenses. Conjugate fault networks record offsets, depressions, or uplifts of individual ore lenses. The acquired understanding of the subsurface from the resulting structural analysis and geologic models improves spatial analysis and can potentially lead to improved mineral exploration targeting. The geologic models and related findings of this contribution serve as a robust foundation for future petrophysical and geophysical studies, and most importantly they can be used as a baseline for integrated multidisciplinary interpretations that are conceptually and quantitatively consistent with all the available data.

This study is part of the Common Earth Modeling (CEM) project financed by Luossavaara Kiirunavaara AB (LKAB), which is thanked for their generous financial support, data sharing, and permissions to publish this study. The authors extend their gratitude to the following individuals from LKAB for their support: Jan-Anders Perdahl, Emma Falksund, Kovács Alpár, Kulcsár Zsolt, Ivan Naumenko, David Monteith, Carlos Santana, Joel Krispinsson, Henrikki Rutanen, and Ulf Bertil Andersson. Special thanks to Laura Lauri and Ian Cope from LKAB for helping during the project, reviewing the manuscript, and providing valuable feedback. We express our sincere appreciation to journal reviewer Laurent Ailleres for providing insightful comments and constructive criticism that greatly improved the manuscript. Additionally, we thank Larry Meinert for his editorial assistance and detailed comments, which contributed to refining the final manuscript. Stefan Luth from the Swedish Geological Survey (SGU) is thanked for his stimulating discussions and valuable inputs, which enriched the study. We acknowledge the use of Leapfrog Geo and extend our thanks to Seequent for their support with the software.

Ervin Veress is a Ph.D. candidate at the Luleå University of Technology (Sweden). His project is dedicated to enhancing the understanding of iron oxide-apatite deposits from deposit to district scale by integrating geologic and geophysical data within an applied geologic modeling framework, utilizing a mineral system approach. He holds an Erasmus Mundus joint M.Sc. degree with a focus on geometallurgy from the University of Liège (Belgium), University of Lorraine (France), and Luleå University of Technology (Sweden).