Implicit Geomodelling of the Merensky and UG2 Reefs of the Bushveld Complex from Open-Source Data: Implications for the Complex’s Structural History

: The Bushveld Complex (BC) is the world’s largest source of platinum group metals. Extensive studies on the complex have focused on its geochemistry, magma and platinum group mineral genesis, mineral characterization and intrusion mechanisms. However, relatively little work has been undertaken on the overall 3D geometry of the complex, which detracts from the adequate contextualization of such studies. Furthermore, the absence of a broader 3D model of the complex does not permit the identiﬁcation of structural trends and mineralization patterns. This contribution details the construction of 3D implicitly-modelled Merensky and UG2 Reefs across the Rustenburg Layered Suite of the BC, using Seequent’s Leapfrog software. Multiple open-source and public-domain data sources and modelling workﬂows were explored to account for disparities in data resolution, data spacing and clustering, and the resolution of model outputs. Key outcomes are (1) a representative, fully-implicit, dynamic geological model of the Merensky and UG2 Reefs over the main chamber of the BC; (2) identiﬁcation of modelled features that augment the current understanding of the BC’s kinematic history and cumulative deformation; and (3) identiﬁcation and analysis of subtle geometrical trends and patterns, such as inter-reef spacing and modelled depths, as well as structural domains that may not have been apparent from numerous, more focused or isolated petrological or geochemical studies. It is anticipated that this baseline 3D model will form the foundation for future, possibly localized, dynamic updates as further information becomes available. The addition of proprietary ( viz. , non-open-source) data, such as 2D seismic sections and 3D seismic surveys, would enhance the overall resolution and quality of such a model and resulting interpretations.


Introduction
The Bushveld Complex (BC) represents the world's largest source of platinum group metals (PGM), which are critical for multiple industrial applications, from automotive catalysts to components in medical equipment. The BC has been the focus of numerous studies, predominantly on whole rock geochemistry and mineralogy, which seek to develop models for PGM genesis and intrusion mechanisms [1][2][3][4][5]. The majority of these studies examine the cumulate layering across various mines and exploration projects, albeit with limited or only very local contextualization within the overall  The UCZ consists of thick anorthosite and noritic anorthosite, with intercalated pyroxenite units and chromite reefs, remarkably regular in mineralogy and thickness, that represent laterally-continuous cyclical units [28]. The UCZ displays well-defined layering, exhibits extraordinary enrichment in PGMs and chromite [15], and represents 50% and 75% of the world's Pd and Pt resource, respectively [26,28]. Given this mineral endowment, concentrated in two main reefs within the CZ, and even though the CZ represents only a relatively small interval of the BC, this interval has been the primary focus of research and mining efforts.
Both the Merensky and the UG2 (Upper Group 2) reefs formed as distinct, persistent layers within the CZ. The Merensky Reef in its "normal" configuration consists predominantly of a coarse-grained, pegmatoidal, feldspathic pyroxenite that contains appreciable phlogopitic mica, varying in thickness from 15 cm to 40 cm, with upper and lower contacts of coarse-grained material typically marked by thin chromitite layers several millimeters thick. Immediately above the top chromitite layer is a brownish pyroxenite known as the Merensky Pyroxenite [1].
The UG2 is a platiniferous chromitite reef that ranges in thickness from 0.4 m to 2.5 m, although occurrences greater than >2.5 m are occasionally observed. This package typically comprises 2 to 9 chromitite layers, intercalated with melanocratic host material that contains a thin chromitite leader [29]. The UG2 footwall comprises either anorthositic units or pegmatoidal pyroxenite [28], with the hangingwall units differing between the EL and WL: The WL hangingwall typically comprises an olivine-enriched pyroxenite, while the EL is predominantly overlain by pyroxenite [29,30]. The position of the UG2 within the cumulate pile is variable, situated between 15 m and 400 m below the Merensky Reef [30][31][32].
Despite its age (2055 Ma [33]), the complex is generally not thought to have undergone any significant metamorphism or deformation other than tectonic subsidence [1,14]. This has been recently challenged by Basson [11], who suggested that deformation of the deeper, central portion of the RLS was responsible for fluid expulsion and generation of the (1) Transvaalide Fold-and-Thrust Event, (2) intrusion of the High-Titanium Suite (HITIS); (3) anticlockwise rotation of the eastern two thirds of the main chamber; (4) thrusting of the northern margin of the main chamber southwards and over the displaced SL and (5) detachment and shearing of the SL from the southern end of the WL.

Geomodelling Theory
Implicit modelling is a spatial interpolation-based modelling methodology. The most widely-used function in 3D geomodelling is the radial basis function (RBF) [34]. Modelled surfaces produced via implicit modelling represent iso-potentials of a scalar field in 3D space [34][35][36]. Generated surfaces are therefore visual representations of isosurfaces based on a 3D scalar field. This allows modelled surfaces to utilize non-contacting data, such as structural trends [34,36]. The advantages of implicit modelling versus explicit modelling are: (1) Modelled surfaces may be pinned locally with hard data, such as drillhole intersections and reef outcrops; (2) Surfaces may be guided by non-contacting data, such as density contrasts from unconstrained geophysical data inversions [37]; (3) Surfaces may be adjusted by user-defined orientation values, in regions where there are no or only sparse available data [38,39]; (4) Extrapolation away from known or hard data, into areas where there is effectively no data, may be achieved intelligently and logically using trends. Such inputs remove, for example, instances where Merensky and UG2 reefs cross one another; and (5) Directional interpretation bias, inherent in section-based explicit modelling, is minimized.
An additional advantage of implicit modelling is the negation of ad-hoc rules during digitization and assembly of lithological solids in explicit, sectional-based digitization, which may be abused in cases where there is complex lithological logging or where the Geologist undertaking the digitization does not or cannot fully understand the overall geological context. Implicit modelling is not, however, without limitations. For instance, multiple, mathematically-induced artefacts have been noted [40][41][42]. These authors noted that local data density has a greater influence on the resultant surfaces, rather than the specific RBF interpolation that is employed. For this contribution, Leapfrog Geo V4.5 was utilized for 3D model construction.

Conceptual Models: Connectivity between Eastern and Western Limbs
The large-scale geometry of the BC has been extensively researched and debated [3,8,35,[43][44][45][46][47][48]. This discussion is constrained by a relatively limited outcrop distribution, a paucity of information at depth-particularly in the central region of the RLS-and the overall scale of the BC. Despite these challenges, several models have been proposed. The EL and WL show marked geochemical, geophysical and textural similarities and, as a result, connectivity between these limbs has been repeatedly proposed (e.g., [1,3,[47][48][49][50][51]). Furthermore, the distribution of intermittent exposures of mafic units, believed to be associated with the BC, suggest the possibility that the complex is also developed or preserved in the center of the two main limbs. As a result, a range of models have been developed for the complex, summarized by Cole et al. [47], each revolving around the degree and/or mode of limb interconnectivity.
Discussions on continuity and geometry towards the central portion of the main chamber of the RLS are hampered by a lack of deep drillholes, forcing a reliance on scattered information from non-PGM/Cr mines and geophysical data inversion (e.g., [11,52]). For instance, the separate dipping sheets model, proposed by Cousins (1959) [51], was initially supported by an apparent lack of RLS distribution at depth and between the two main limbs. The same study interpreted intrusions with separate vertical feeders (see also Viljoen, 1999 [52]), such as those that emanate from the edges of a mantle diapir penetrating the crust. Older gravimetric studies and geoelectric surveys [7,[53][54][55] surmised a lack of mafic units between the Dennilton Dome and the Rooiberg Fragment ( Figure 1). Additional geophysical work by Cheney and Twist (1991) [54] appeared to support the absence of RLS at depth, between the EL and WL. Models of individual, inward-dipping sheets and discrete or separate limbs have also been presented by several authors [7,[55][56][57][58].
The intrusion of a relatively flat sill into a broadly folded basin was proposed by Gruenewaldt, (1979) [57], while Eales et al., (1993) [26] surmised that these sheets were initially emplaced horizontally. In the main chamber, the RLS either shallowly cross-cuts the upper Pretoria Group or follows its contact, while the later Lebowa and Rashoop Suite granites intruded into the upper portions of the RLS and Rooiberg felsites [10]. Several authors have suggested that deformation preceded the complex's intrusion and that strain or folding may be attributed to the subsidence of the encompassing sedimentary basin [51,59,60]. Regional geological mapping outlined large fold structures and large country rock xenoliths in the floor of the complex [61][62][63]. Tectonic fabrics and signs of dynamic metamorphism, including bedding-parallel fabrics in Pretoria Group pelites, were recorded. However, Button (1973) [64], Blain (1975) [65], and Walraven (1981) [66] noted that Transvaal Supergroup units within the vicinity of the complex appear to be "relatively undeformed and display the original basin shape rather well". More recent studies support inter-limb continuity-or at least greater inter-limb continuity-at depth [1,3,9,59,67,68]. Geophysical data modelling by Cole et al. (2013) [47] indicated a significant volume of BC material between the EL and WL, west of the Wonderkop and Stofpoort Faults (fault positions indicated in Figure 2), which supports the concept of continuous mafic units at depth. Cole et al. (2013) [47] also noted that magnetic anomalies of the RLS are heavily influenced by magnetite-rich UCZ units, which cause the majority of discernible positive magnetic anomalies and identifiable layering.  Work by Kinnaird and McDonald (2005) [61] and Kinnaird et al. (2017) [69] indicated that the NL is likely detached from the main part of the main chamber by the ENE-trending Thabazimbi Murchison Lineament, which may have comprised a dyke-like feeder for Bushveld and other magmas, resulting in an impediment between the NL and the main body of the complex. Gravity models by Finn et al. (2015) [48] indicate moderately thicker and/or high(er)-density materials at fault intersections, within the western and central TML, which have been interpreted as feeders, although this interpretation may rely on the assumption that deformation in these areas was fairly limited.
There are four main geometrical models for the main chamber: (1) The EL and WL are separate or discrete entities [51]; (2) The complex is continuous at depth, viz., the connected lobes model as suggested by [1,3,9,70]; (3) The complex is mostly continuous, with local disruption by major structures or domes; the separate dipping sheets model [7,55] or (4) The central part of the RLS has been extensively displaced by southwards-directed thrusting, the EL and WL at depth have undergone a degree of thrusting and folding and the eastern portion of the main chamber has undergone anticlockwise rotation of approximately 35 • [11]. Basson (2019) [11] proposes that the geometry in the central portion of the main chamber is structurally complex and does not reflect the original emplacement geometry.

Materials and Methods
A broad overview of the modelling methodology will be addressed, outlining the workflow and development of the BC implicit geomodel from data incorporation to results and analysis. In the case of 3D geomodelling in Leapfrog Geo (Software version 4.5.1, Seequent, New Zeeland), this procedure follows a standardized, basic workflow: (1) Establishment of a volume of interest; (2) Incorporation of data and its quality control; In Leapfrog Geo, geological relationships are established as a series of conditional assembly rules that are geologically-phrased and must be geologically logical. An example of this is the establishment of a chronological sequence dictating how the stratigraphic or tectonostratigraphic sequence should be assembled. However, this methodology is not without its limitations, and the construction of surfaces and volumes produced during syn-deformational emplacement is difficult without user input, particularly in areas of low data density. Once contacting and age relationships between geological units have been defined, geological parameters are established for each surface. Surfaces may be modified according to additional data or observed deformation style(s). Structural trends may be incorporated via a number of surface construction tools, allowing for greater surface control and the generation of geologically-realistic surface geometries. Once contact surfaces have been used to subdivide the 3D volume into subdomains, resultant outputs are amenable to iterative interrogation, adjustment of parameters and subsequent modification by the insertion of additional controls. Re-running the model, using these updated surfaces, automatically produces updated 3D volumes or wireframes.

Model Construction
The volume of interest ( Figure 2) encompasses the bulk of public-domain data across the RLS of the main chamber of the BC, with its eastern, western and southern margins defined by a buffer region of 500 m outwards from surface UG2 and Merensky Reef outcrops. The northern extent is a 500 m buffer, northwards of structures interpreted to form part of the TML [11,71]. The TML is a well-documented, regional tectonic boundary which, along with the Mohlapitsi Fold-and-Thrust Belt, the Welgevonden/Southern Fault System and the Zebedelia Fault, trends ENE-WSW for approximately 600 km across the Kaapvaal Craton, thereby representing a broad, regional fault-and-shear zone that was active syn-to post-BC emplacement [11,54,72]. Using these constraints and a modelled depth of −4500 MASL, the model has a lateral extent of approximately 400 km by 330 km and a total volume of approximately 270,000 km 3 . The central and northern areas of the main chamber are largely overlain by sediments of the Waterberg and Karoo Supergroups, with few identifiable BC outcrops or occurrences, although several deep drillholes have indicated the presence of RLS units [73]. The geological model presented herein is constructed in the WGS 1984 Universal Transverse Mercator (UTM) 35S coordinate system, as the majority of available data falls within this UTM zone.

Structural Network Incorporation
A core component of 3D geomodelling is the early establishment of a representative fault network [38,70], which ideally exhibits a consistent internal deformation history. Key faults are based on the structural framework derived from Basson (2019) [11]. This compiled fault network and additional surface and underground mapping structural data were incorporated into 3D space for fault construction, which further sub-domains field mapping, geophysical, and stratigraphic/drillhole datasets. A total of 63 fault traces were incorporated into a cohesive fault network, re-constructed in Leapfrog Geo. Fault characteristics, such as dip and dip orientations, were captured from available literature and explicitly characterized in terms of their chronology or crosscutting relationships. Thirteen regional faults were selected to sub-domain the 3D volume into 23 individual sub-volumes ( Figure 3).
Minerals 2020, 10, x FOR PEER REVIEW 9 of 24 crosscutting relationships. Thirteen regional faults were selected to sub-domain the 3D volume into 23 individual sub-volumes ( Figure 3). A total of 63 key regional structures were identified and modelled in 3D. Of these, 13 were activated within the 3D modelling software (Leapfrog), to create 23 discrete fault-bounded domains. Merensky Reef and/or UG2 Reef surfaces were not modelled in several fault domains due to a paucity or lack of outcrop or downhole data. Static structural trends were utilized for initial reference surface construction, which varied across fault domains. These structural trends were derived from stereonet analysis and contoured averages of cumulate layering measurements and, if available, regional bedding orientations (stereonets are equal area, lower hemisphere). Data sourced and modified/augmented after Basson (2019) [11].

Surface Construction Technique
Limitations in software capability necessitated a three-stage process to generate the key marker reefs. The modelling workflow is detailed in Figure 4: (1) Construction of reference surfaces for both Merensky and UG2 reefs from mapping data, drillhole-derived downhole reef intersections and surface reef contours; (2) Construction of structural trends, incorporating trends and surfaces produced in unconstrained geophysical data inversions (see Basson, 2019 [11] for a description of these); and (3) Construction of final surfaces. A total of 63 key regional structures were identified and modelled in 3D. Of these, 13 were activated within the 3D modelling software (Leapfrog), to create 23 discrete fault-bounded domains. Merensky Reef and/or UG2 Reef surfaces were not modelled in several fault domains due to a paucity or lack of outcrop or downhole data. Static structural trends were utilized for initial reference surface construction, which varied across fault domains. These structural trends were derived from stereonet analysis and contoured averages of cumulate layering measurements and, if available, regional bedding orientations (stereonets are equal area, lower hemisphere). Data sourced and modified/augmented after Basson (2019) [11].

Surface Construction Technique
Limitations in software capability necessitated a three-stage process to generate the key marker reefs. The modelling workflow is detailed in Figure 4: (1) Construction of reference surfaces for both Merensky and UG2 reefs from mapping data, drillhole-derived downhole reef intersections and surface reef contours; (2) Construction of structural trends, incorporating trends and surfaces produced in unconstrained geophysical data inversions (see Basson, 2019 [11] for a description of these); and (3) Construction of final surfaces.

Reference Surface Construction
Reference surfaces are required for the application of software-related tools on such widelydistributed and variably-scaled datasets. Foremost amongst these tools is the ability of the software to utilize non-contacting structural data to inform surface construction, using Leapfrog's Structural Trend tool. This tool produces a scalar field where the gradient of that field is influenced by structural data. Reference surfaces for the UG2 and Merensky Reefs were constructed, utilizing contacts derived from drillhole logs, surface mapping and reef contours. These reference surfaces utilize static guiding trends (fault domain-specific trends displayed in Figure 3), in turn derived from various technical reports, maps and academic studies. Given the data distribution, the deposit surface tool allows for construction of the most robust reference surfaces with a minimum of explicit controls, for the UG2 and Merensky reefs. The deposit surface accomplishes this by explicitly using specified contacts or intervals, even over great distances.
Control points were then derived from the modelled UG2 and Merensky Reef reference surfaces, which denote both the position and orientation of the related surface triangulation [74]. Given the scale of the model, these control points required simplification or declustering prior to downstream use. Declustering is the procedure whereby large datasets are reviewed, sets of typical orientations identified and averages, per set, calculated over specified distances or volumes. The Micromine software package (Software version 2020.5, Micromine, Nedlands, WA, Australia) was utilized to decluster the resultant reference surface control point files within a 1 km buffer zone, as these tools are not available within the Leapfrog modelling package. This provided a calculated mean orientation at a centroid for all points within a 1 km radius. These control points were subsequently reimported into the Leapfrog modelling environment.

Reference Surface Construction
Reference surfaces are required for the application of software-related tools on such widely-distributed and variably-scaled datasets. Foremost amongst these tools is the ability of the software to utilize non-contacting structural data to inform surface construction, using Leapfrog's Structural Trend tool. This tool produces a scalar field where the gradient of that field is influenced by structural data. Reference surfaces for the UG2 and Merensky Reefs were constructed, utilizing contacts derived from drillhole logs, surface mapping and reef contours. These reference surfaces utilize static guiding trends (fault domain-specific trends displayed in Figure 3), in turn derived from various technical reports, maps and academic studies. Given the data distribution, the deposit surface tool allows for construction of the most robust reference surfaces with a minimum of explicit controls, for the UG2 and Merensky reefs. The deposit surface accomplishes this by explicitly using specified contacts or intervals, even over great distances.
Control points were then derived from the modelled UG2 and Merensky Reef reference surfaces, which denote both the position and orientation of the related surface triangulation [74]. Given the scale of the model, these control points required simplification or declustering prior to downstream use. Declustering is the procedure whereby large datasets are reviewed, sets of typical orientations identified and averages, per set, calculated over specified distances or volumes. The Micromine software package (Software version 2020.5, Micromine, Nedlands, WA, Australia) was utilized to decluster the resultant reference surface control point files within a 1 km buffer zone, as these tools are not available within the Leapfrog modelling package. This provided a calculated mean orientation at a centroid for all points within a 1 km radius. These control points were subsequently reimported into the Leapfrog modelling environment.

Construction of Structural Trends
The Structural Trend tool generates an ellipsoidal anisotropy that varies in direction and strength according to an input triangulation, surface, or structural measurement. These ellipsoid ratios determine the local weighting of data informing surface construction. These structural trends may be viewed or envisaged as a series of "flattened spheres" in 3D space, where the orientation of a sphere provides the direction of the anisotropy, while the size of a sphere or disk is proportional to the local anisotropic strength. A structural trend was constructed per fault domain, to inform final surface construction trends within each individual domains. These structural trends incorporated field mapping, drillhole contacts, reef contours and geophysical inversion products [37,52,74], which all have disparate or different scales, densities and coverage [75].
As indicated previously, contoured isoshells of gravity contrasts, representing inferred depths and geometries of key contacts of the RLS [11], were used to guide RBF-generated UG2 and Merensky surfaces at depth, away from borehole-derived contacts, reef contours or reef traces, similar to the statistical approach used by Gupta and Sutcliffe (1990) [76]. Incorporating the gravity inversion products required additional upfront processing. Unconstrained gravity inversion contours were used to construct a series of reference surfaces, from which control points, which indicate the spatial position and orientation of key contacts, could be extracted. Multiple surface iterations were reviewed to ensure that the most geologically-representative surface was selected for downstream modelling.
Gravity highs were identified where the RLS approaches the surface towards the West in the WL and towards the East in the EL. Upon review of the control point file, it was noted that several control points indicated a down-dipping geometry just beyond RLS surface outcrops to the West and East of their respective limbs. This apparent down-dipping geometry was due to the decrease in the density signal beyond the surface exposure of RLS and consequently these non-representative down-dip control points were selectively removed.
To ensure that both unconstrained inversion data as well as regional field mapping trends were effectively utilized, two regional form interpolants were generated. A form interpolant is a representative surface, or series of surfaces, that emulates broad structural anisotropy from strike/dip trends in 3D, which may be used as an input for structural trend generation [74,77,78]. Two separate form interpolants were generated: One derived from the unconstrained gravity inversion and the other from compiled field mapping data. Unconstrained gravity inversion data required post-processing, where resultant isoshells were converted to a series of polylines and orientation points in 3D space. Once generated, these 3D elements were then used to define their resultant form interpolant. The resultant surfaces for these form interpolants were used as primary inputs for fault domain-specific structural trends, resulting in a "blended" structural trend that incorporates all available structural data ( Figure 4).

Construction of Final Surfaces
Final UG2 and Merensky contact surfaces were then generated (Figures 4 and 5), incorporating declustered control points and initial hard lithological contacts/surface traces, with their geometries dictated by fault domain-specific structural trends. These surfaces were constructed initially as intrusion contact surfaces, which allowed for the implementation of the structural trend tool.

Results
Modelled contact surfaces, including contoured reef depths, surface dip values and inter-reef separation distances, may have implications for the degree of erosion of reef footwalls. Reef depth maps for the Merensky Reef ( Figure 6a) and UG2 Reef (Figure 6b) were contoured at a vertical spacing or interval of 10 m. WL Merensky and UG2 Reef depth contours indicate substantial increases in depth (−3000 m MASL) towards the center of the complex, east of the Crocodile River Fault, although one needs to take into consideration the effects of extrapolation into areas with sparse data constraints, where surface construction defers to the inversion of regional geophysical data. This substantial increase in depth coincides with areas overlain by the Lebowa Granite Suite. Modelled Merensky and UG2 reefs shallow to the east of the Crocodile River Fault, indicating a trough-like depression to the west, possibly accommodated by west-up vertical displacement along the Crocodile River Fault. The Merensky Reef also shows significant variations in depth across the Brits Graben Fault, although this structure appears to have affected the depth and/or depth variation of modelled UG2 in this area to a lesser extent.

Results
Modelled contact surfaces, including contoured reef depths, surface dip values and inter-reef separation distances, may have implications for the degree of erosion of reef footwalls. Reef depth maps for the Merensky Reef ( Figure 6a) and UG2 Reef (Figure 6b) were contoured at a vertical spacing or interval of 10 m. WL Merensky and UG2 Reef depth contours indicate substantial increases in depth (−3000 m MASL) towards the center of the complex, east of the Crocodile River Fault, although one needs to take into consideration the effects of extrapolation into areas with sparse data constraints, where surface construction defers to the inversion of regional geophysical data. This substantial increase in depth coincides with areas overlain by the Lebowa Granite Suite. Modelled Merensky and UG2 reefs shallow to the east of the Crocodile River Fault, indicating a trough-like depression to the west, possibly accommodated by west-up vertical displacement along the Crocodile River Fault. The Merensky Reef also shows significant variations in depth across the Brits Graben Fault, although this structure appears to have affected the depth and/or depth variation of modelled UG2 in this area to a lesser extent.
UG2 and Merensky Reefs increase in depth (−1700 m MASL) towards the Waterberg Basin, north of Bela-Bela, in an area dominated by faults and thrusts south of the TML, although this increase in depth is more apparent in modelled UG2, possibly suggesting thrust-related splaying. Modelled Merensky and UG2 reefs shallow significantly towards the Marble Hall Fragment, Dennilton Dome and the Steelpoort Fault to the east. Modelled reefs also appear to be locally disrupted by minor domes within the Eastern Limb (e.g., Signal Hill/Boschpoort/Makgane Domes), although reef depths appear to be relatively consistent E and NE of the Burgersfort Bulge. Reefs shallow significantly to the SE of the Steelpoort Fault, north of Booysendal Mine.
Modelled inter-reef separation is presented in Figure 7. Zones of anomalously high inter-reef separation were identified (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12). Upon review, it was established that several of these anomalous inter-reef separation highs (1-6) may be modelling artefacts, caused by disparities in data availability for the modelled Merensky Reef and the possible effects of over-extrapolation. One of these artefacts is an artificially-uplifted or higher Merensky Reef surface, due to the modelling interpolant attempting to connect data that are positioned at a greater average height than what is appropriate for a given fault-domain across that modelling domain.    Modelled inter-reef separation is consistent along the majority of the WL, particularly in proximity to reef outcrop, with a typical separation of between 10 m and 50 m. The northern part of this limb is known to host many potholes (see Smith and Basson, 2006 [79], circular or ellipsoidal depressions that may range from several meters to several hundred meters in diameter [11,80,81], although at the scale of the model presented herein, inter-reef separation is partially smoothed out. This pattern is locally disrupted east of the Elandsdrift Fault, towards the Krokodildrift Section, where inter-reef separation dramatically increases and in some areas is >400 m (e.g., Point 7). Separation apparently increases around the boundaries of the Pilanesberg Complex, ranging between 100 m and 200 m (Point 8). The model suggests significant inter-reef separation (400 m to 500 m) between the Middellaagte Graben Bounding Fault and the northern portion of the Crocodile River Fault (Point 9). West of the Dennilton Dome (Point 10), to the south of a regional dyke, modelled inter-reef separation increases and may be as high as 350 m. Similar to the WL, the modelled EL Merensky and UG2 Reef separation is largely consistent in proximity to reef outcrop, with separation distances constrained to between 10 m and 100 m. Reef separation increases significantly to the northeast of the Dwarsrivier Splay, in close proximity to Kennedy's Vale Mine (Point 11), with modelled inter-reef separation ranging from 300 m to 500 m. Similarly, anomalously high modelled inter-reef separation is again evident south of the Laersdrift Fault (Point 12). Moderate inter-reef spacing is apparent in the vicinity of the Marble Hall Fragment (100 m to 200 m), increasing dramatically towards the SE, although high inter-reef separation appears to dissipate or is inconsistent near several regional dykes (Point 13), which may represent areas of fault-related drag.
Modelling artefacts, particularly those relating to modelling at depth, are to be anticipated in a geological model of this scale and complexity [79], using widely separated, locally disparate data sets at a variety of scales. As the spatial RBF function attempts to generate a smooth surface through control points, areas of sparse data are likely to produce over-simplified surfaces or surfaces that are locally unrepresentative of geological features or processes, particularly where the tectonostratigraphy changes its geometry abruptly and/or where data are sparse. Furthermore, implicitly modelled surfaces and resultant solids may extrapolate unrealistically outwards and away from their constraining data, towards model boundaries. These instances have to be noted and mitigated, where possible. Notwithstanding this, implicit 3D geomodelling over such a large area is eminently possible, with the application of a suite of modelling assumptions and controls.
Dip maps of the modelled Merensky and UG2 Reefs are displayed in Figures 8 and 9. Areas of dip maps, in the context of the Bushveld Complex, which show steep dips may denote slumping or potholing of reefs. Dips steepen towards surface reef outcrops across both the WL and EL of the BC. Anomalously steeply-dipping modelled Merensky and UG2 Reef occurs towards the northern and southern terminations of the WL (see cross-section on Figure 9). These steep dips possibly extend eastwards from the Crocodile River Fault system, for approximately 5 to 6 km (Point 1). Extensive or pervasive steepening of the modelled Merensky and UG2 Reefs is observed to the north of the Welgevonden Fault Zone, although such steep dips are more apparent in the contoured UG2 Reef dip map (Point 2). Anomalously steep dips are observed to the north of the Wonderkop Fault splay and extend southwards towards the Stavoren/Makekaan Fragment (Points 3). Broad zones of almost horizontal reef are noted towards the apparent or inferred center of the RLS in the main chamber (Point 4), and towards the southern portion of the EL (Point 5), forming broadly concentric "domes" of low average dips, although this may partially be due to the lack of data in these areas. A general steepening of the modelled UG2 Reef surface and cumulate layering (<50 • ; stereonets in Figure 5) is observed towards the Zebedelia Fault and the northern edge of the EL (Point 6).

Discussion
Multiple data sources, modelling workflows and solutions were explored during this study to account for the disparities in data resolution, data spacing/clustering, model scale and output. Resultant geometries and interpretations are nevertheless dependent on data distribution and quality, which is the case for any 3D model. Added to this is the modeler's explicit controls, choice of modelling technique, and available tools within software packages such as Leapfrog and Micromine. However, the position and incorporation of additional data will have a significant impact on modelling results and such data, be it "softer" geophysical data or "hard" drillhole data, are likely to have much greater influence towards the center of this geological model where currently-available data are sparse. Seismic sections, which may have already been obtained by mining and exploration companies, would be particularly helpful in refining the 3D model, which, together with the workflow, was constructed with the intent of iteratively incorporating future additional data, in parallel with the activation of existing or even additional faults as required. This, in conjunction with more local-or mine-scale comparisons, would provide a series of increasingly better-constrained outputs that may be similarly subjected to interpretation, benefitting, in particular, where higher-resolution data have an impact on smaller fault domains wherein information is currently relatively sparse or diluted.
The northern portion of the EL appears to be absent to the west of Messina Mine. The modelled Merensky and UG2 Reefs shallow, and then steepen north-upwards, towards and along the TML-Zebediela-Welgevonden/Southern Fault systems. In addition, the EL becomes progressively more faulted, evidenced or accompanied by rapid shallowing of both reefs (Figures 8 and 9), towards the northern margin of the main chamber and in proximity to the left-lateral Stofpoort/Wonderkop Fault system. This apparent uplift and steepness along the northern model margin, in conjunction with the development of the Stofpoort/Wonderkop Fault System(s), likely developed in response to N-S compression of the complex [11,72]. This southward vergence and asymmetry are also evident from the Droogekloof Thrust-Southern Boundary Fault system near the Thabazimbi/Rooiberg area. The variation in modelled inter-reef separation in this area may represent the effects of a forward-propagating thrust system that acted through the central part of the main chamber, as proposed by Basson (2019) [11].
Elevation plots suggest deep, relatively undeformed, relatively flat-lying Merensky and UG2 Reef geometries east of the Crocodile River Fault, with Merensky and UG2 Reefs apparently shallowing eastwards. Notwithstanding this, the RLS to the west of the Crocodile River Fault appears to have accommodated N-S flattening, if one assumes an initially more rounded, flatter-dipping magma chamber, by developing steeper northern and southern edges (see cross-section in Figure 9). The southern side of the section on Figure 9 suggests that the presence of the Johannesburg Dome, as proposed by Basson (2019) [11], was instrumental in steepening reef dips. Anomalous inter-reef separation, significant depths and steep bedding angles at the northern and southern terminations of the WL likely indicate that southward vergence due to N-S compression is more clearly or unambiguously evident in these parts of the WL, and particularly in proximity to the termination of the Krokodildrift Section. Low-angle bedding or cumulate-layer parallel shears, recorded in underground platinum mines on the southern part of the WL [11,21,79,80] and sheared contacts between tectonically-interlayered Transvaal Supergroup metasediments and Bushveld Complex-related Norite [11,49,50,67] only occur to the east of the section shown on Figure 9. In contrast, the section shown as an inset to Figure 8, across the EL, is markedly different from that across the WL. This shows significant steepening in the north, but an uplifted area in the vicinity of the Phepane Dome, and relatively flat-lying geometries across the remainder of the section to the south. This suggests that N-S strain was accommodated in a different manner in the EL, combined with the rotation, proposed by Basson (2019) [11] of the eastern 2/ 3rds of the main chamber.
These observations, in conjunction with the homogenous reef dip in the area, indicate that the Crocodile River Fault has likely acted as a regional strain barrier, effectively partitioning strain between the WL, which shows that N-S flattening was accommodated by the upliftment and steepening of its northern edge, and impingement and steepening, against the Johannesburg Dome, of its southern edge. The remaining 2/3 rds of the main chamber accommodated N-S flattening, by anticlockwise rotation and the probable development of a forward-propagating thrust splay system. As a result, the WL also appears to be deeper, i.e., it has been relatively more down-faulted along the Crocodile River Fault. While still exposed at surface it is more "bowed" compared to the EL, which accommodated the N-S imposed strain by rotating.
Another pattern emerges from the inter-reef separation map: there is a zone of apparently low inter-reef separation, approximately 50 km SW of the Dennilton Dome. Deformed Transvaal Supergroup inliers occur in the BC, the most prominent of which are the Crocodile River Dome (western Transvaal Basin), the Dennilton-Marble Hall Dome and Stavoren Fragment (eastern Transvaal basin), which consist predominantly of deformed lower Transvaal units that underwent low grade metamorphism [81,82]. Positions and distributions of low dip angles, in conjunction with the positions of various domes towards and in the EL, may again support the concept of a forward propagating thrust system to the east of the Crocodile River Fault, which could have transported fragments (including "domes") of the northern margin of the main chamber southwards.
In the area between the southern terminations of the EL and WL, moderate to low reef depths, moderate to low reef dips and low inter-reef separations emerge from the 3D modelling exercise, suggesting that reefs from the southern margin of the main chamber, in this area, are likely to have experienced layer-parallel shearing during southerly vergence or have been transposed into parallelism with thrusts, probably within a forward-propagating thrust system that stems from or is rooted in the TML. In addition, Basson (2019) [11] proposed that the offset between the SL and the southern termination or tip of the WL was right-lateral, accommodating an apparent 35 • anticlockwise rotation of the volume to the east of the Crocodile River Fault, with the initial BC body behaving in a ductile manner at depth. Strain accommodation near-surface was brittle or brittle-ductile, in a collection of mutually cross-cutting "relay faults" in and around Cullinan Mine. The coincidence of deep reef and the overlying presence of the Lebowa Granite Suite is apparent from this study.

Conclusions
Modelled reef geometries support a model of various modes of stress accommodation, and thereby modes of strain, within the main chamber of the BC to the south of the TML. The WL shows steepening at its northern and southern extremities, against the TML and Johannesburg Dome, respectively, resulting in downward "bowing" of reef surfaces and E-up displacement along the Crocodile River Fault and related, parallel right-lateral faults in the vicinity.
The EL, on the other hand, probably preserves the original geometry of the intrusion somewhat better, suggesting very flat dips of the original cumulate pile. However, the EL is not undeformed. It displays steepening of its northern margin, against the TML, but it appears that strain was accommodated by a degree of (real or partially apparent) left-lateral shearing along the TML and anticlockwise rotation of the eastern 2/3 rds of the main chamber [11]. The most likely structural configuration of the area immediately to the east of the Crocodile River Fault, which would have aided in the preservation of flat reef dips, is a southward-and forward-propagating thrust system, evidenced by entrained BC material within the Cullinan Mine [11,73]. The Crocodile River Fault therefore acts as a strain delimiter, separating modes or styles of deformation.  Acknowledgments: The first author thanks Ian Basson for supplying the substantial underlying database, as well as his time and expertise during the completion of this study, without which this project would not have been possible. The author acknowledges the financial, logistical and software support provided by Tect Geological Consulting. The assistance of Marietjie Schalekamp from the South African Council for Geoscience (CGS), in the sourcing of drillhole data, is further acknowledged.