Constraining 3D geometric gravity inversion with 2D reflection seismic profile using a generalized level-set approach: application to Eastern Yilgarn craton

One of the main tasks in 3D geological modelling is the boundary parametrization of the subsurface from geological observations and geophysical inversions. Several approaches have been developed for geometric inversion and joint inversion 15 of geophysical datasets. However, the robust, quantitative integration of models and datasets with different spatial coverage, resolution, and levels of sparsity remains challenging. One promising approach for recovering the boundary of the geological units is the utilization of a level-set inversion method with potential field data. We focus on constraining 3D geometric gravity inversion with sparse lower-uncertainty information from a 2D seismic section. We use a level-set approach to recover the geometry of geological bodies using two synthetic examples and data from the 20 geologically complex Yamarna terrane (Yilgarn craton, Western Australia). In this study, a 2D seismic section has been used for constraining the location of rock unit boundaries being solved during the 3D gravity geometric inversion. The proposed work is the first we know of that automates the process of adding spatially distributed constraints to the 3D level-set inversion. In many hard-rock geoscientific investigations, seismic data is sparse and our results indicate that unit boundaries from gravity inversion can be much better constrained with seismic information even though they are sparsely distributed within the model. 25 Thus, we conclude that it has the potential to bring the state of the art a step further towards building a 3D geological model incorporating several sources of information in similar regions of investigation.

only allows us to define an arbitrary number of density contrast units free of shape limitation but also allows us to add sparsely distributed low uncertainty data to constrain the gravity inversion. We use information from the seismic section as low 60 uncertainty data as in large-scale studies borehole data might not necessarily be available for constraining purposes. For more details about the methodology, we refer readers to (Li et al., 2017;Raponi et al., 2017;Tai & Chan, 2004;.
In this study, we focus on constraining surface gravity data that possesses good lateral resolution with a 2D reflection seismic profile that traverses the study area. In this manner, sparsely distributed low-uncertainty data from the seismic section can be utilized to guide 3D gravity inversion results. The proposed work is the first we know of that quantitatively integrates the 65 geometries from the seismic section into 3D gravity inversion. Our results suggest that the proposed approach has the potential to bring the state of the art a step further towards automatically building a 3D geological model which is consistent with available geological and geophysical datasets.
The remainder of this paper is organized as follows. We first provide a summary of the generalized level-set method we use and show its applicability to gravity inversion constrained by sparser seismic data and geological knowledge. We then apply 70 the proposed approach on two 3D synthetic datasets to show the proof-of-concept. Subsequently, we present a case study at the geologically complex Yamarna terrane in Yilgarn Craton, Western Australia. We use the 3D surface gravity datasets and a 2D seismic profile available in this area to apply the constrained inversion approach. Finally, the performance of the constrained level-set approach on the synthetic and field studies and the resulting models are reviewed and assessed in the discussion section. 75 2 Method

Generalized level-set method
We use the generalized level-set inversion formulation introduced by Giraud et al (2021a) for our constrained inversion problem. We extend their work to the use of sparse seismic constraints and topological rules. A summary of the method is provided below. We refer to Giraud et al (2021a) for details about the mechanics of the inverse problem and use the same 80 notation as these authors in what follows.
First, a geologically plausible starting model is defined with density contrasts ∆ =1,.., assigned to each of the N rock units.
Then for each unit, signed-distance values to interfaces ( ) are calculated using the fast marching method (Sethian, 1999), where the outline of a given unit is defined by = 0. A smeared-out Heaviside function is then defined to transform the signed-distance functions to a multinary structure. This allows us to generate a physical property model (  (1). In this framework, residuals of the calculated and observed datasets ( = [ − ] ) are minimized during the inversion. The data misfit term to optimize can be written as follows: 95 Where Ψ r represents the misfit function.
The misfit function is minimized iteratively to solve for changes in the signed-distance function ( ). To stabilize the inversion problem and incorporate prior information, regularization terms are added to Eq. (1) as discussed below. We use an updated form of the regularizations as spatial constraints to the inversion problem based on seismic information using the 100 aforesaid method. To encourage geological plausibility, we then introduce the utilization of topological rules on the resulting model during inversion. The utilization of seismic information to define regularizations is introduced below.
In such cases, the problem is usually regularized by favouring structures with the shortest interface overall length (2D case) or 105 smallest surface (3D case). Smoothing the geometries characterized by the zero level-set functions by minimizing the length, generates shapes with the smallest area and regularizing these inversion problems can be limited to the specific shape of units that can introduce a bias towards unrealistically simple geometries. Using known density contrast values for the model parameterization in these approaches significantly reduces the non-uniqueness of the inversion. In addition, regularizing the inverted model using prior information may reduce the non-uniqueness of the inverse problem further. In the level-set 110 inversion scheme we use, prior information can be appended as regularization terms (Ψ t ) in the same fashion as that smallness terms regularize inversion problems (Calvetti et al., 2000). This term is appended to the misfit function as below: Where; Where Ψ t tends to minimize the total update in the signed-distance ( ) and acts as a smallness term that encourages the signed-distance update ( ) of rock units to reach specific values stored in = ( , 1 , 2 , … , ). Here, is a global regularization 1 × MN vector, and is a local regularization N × MN matrix that has =1,2,… as rows. Both and 120 encapsulate prior information for the inversion.
Inverting for the total changes in the signed-distance value at each model update and controlling these changes through regularizations enables us to extend the approach for the constrained inversion problems. In the next subsection, we introduce the method we developed to incorporate 2D information into the regularization scheme presented here, using the case of 2D seismic data. We show that updating regularization terms based on low-uncertainty datasets extends the application of global 125 and local regularizations to act as global and local constraints.

Translating seismic information to spatial constraints
Global ( ) and local ( =1,…, ) regularization terms are appended to the sensitivity matrix of the level-set method as in Eq. (5) as uniform vectors in a system of equations solved in the least-squares sense. They tend to stabilize changes of signeddistance functions as damping terms. This supports the capability of regularization terms to perform the constraining process 130 using low-uncertainty datasets such as seismic as weighing terms. We propose that adding unevenly distributed weights to the regularizations based on information from seismic in depth can act as constraints for the level-set inversion. Therefore, as long as the interpreted or inverted sections from the seismic are included in vectors with the same size as the model, any primary information from pre-existing modelling such as seismic sections can be translated to weighting terms to define constraints.
Strategies for transferring the knowledge from a seismic section to a weighting term that can be used as sparse constraints can vary depending on the seismic data and availability of other datasets. Overall, the interpretation from a vertically extended seismic profile can be assumed as an interpreted sample geological model. In Fig. 1, we assume that interpretation has led to https://doi.org/10.5194/se-2021-65 Preprint. Discussion started: 16 June 2021 c Author(s) 2021. CC BY 4.0 License. a conceptual geological model comprising of N lithological units. Therefore, 2D seismic interpretations are treated in the same fashion as a 2D geological surface map to constrain the gravity inversion at depth. Constructing the weighted regularization 140 terms or constraint terms based on available interpretation is necessary before starting the inversion. All arrays of the model that lie within the sample section are weighted accordingly for all lithologies as =1,…, (all of dimensions 1 × M) which then construct the regularization terms as in Eq. (6). Thus, the values of and =1,…, are adjusted locally accordingly with seismic interpretations to favour the preservation of interpreted units along the seismic section. The structure of weighting vectors from seismic interpretations within regularization terms is as follow: 145 Where, is a zero vector of dimensions 1 × M. Figure 1 shows the weighting matrix for a given lithology (Fig. 1a) and its location within regularization terms. In this case, in a matrix of one value with the same size as the model, the entire extracted section along the seismic profile is weighted accordingly for all rock units. These spatially distributed weighting sections (matrices) are then transferred to vectors and 150 applied as global and local constraints to the previously introduced regularization terms. It is important to increase the values of weights to reduce model updates around the seismic interpretations and to suppress the changes of the interface locations in low-uncertainty areas.
As the calculation of the sensitivity matrix ( ) is limited to the vicinity of interfaces defined by , the search space for applying the constraints can also be limited to neighbouring cells around interfaces. This indicates that regularizations can also 155 be adjusted over boundaries between units rather than the whole lithological units. This can be advantageous in cases where extracting the array locations along the boundaries is straightforward. Detected unit boundaries from other sources of information (geological and geophysical) can then be directly transferred to weighting matrices to be used in the constrained level-set inversion. Boundaries' parametrization from interpretations or inversion (such as acoustic impedance inversion) of seismic data after depth conversion and projection onto the mesh used for the inversion can be directly used as regularizations. 160 It is then plausible if interpreted boundaries from regional seismic studies or from inversion are digitized as weighting matrices

165
Distribution of constraint matrix of from lithology 2 for the entire unit (a) and for the corresponding boundaries (b). is shown only for the corresponding 2D section.
As the application of sparse constraints enforces local restrictions to the evolution of signed-distance values, the resulting model from ensembles of signed-distance functions for different lithologies also requires constraints to enforce small-scale topological rules accordingly with geological knowledge and to ensure that certain configurations do not occur during 170 inversion. This is covered in the next subsection.

Enforcing topological rules
Topological constraints play a paramount role in geo-modelling processes (Burns, 1988;Pellerin et al., 2017;Perrin and Rainaud, 2013;Thiele et al., 2016). They refer to properties of a model that are changing during the perturbation of a model.
Topological relationships can be defined over discontinuity networks (fractures or faults) (Sanderson and Nixon, 2015), 175 lithological units, or unconformities (Jessell et al., 2010) and are decisive components in 3D modeling process. In recent years, these relationships are computationally being considered either explicitly or implicitly (Pakyuz-Charrier et al., 2018a;Pakyuz-Charrier et al., 2018b) and as constraints for probabilistic 3D models (De La Varga et al., 2019) in the context of 3D geological modelling. Furthermore, in the geophysical inversion context, they have been sometimes addressed by means of uncertainty https://doi.org/10.5194/se-2021-65 Preprint. Discussion started: 16 June 2021 c Author(s) 2021. CC BY 4.0 License.
quantifications during inversion and joint-inversion problems (Giraud et al., 2019b(Giraud et al., , 2017Wellmann et al., 2018) or as post-180 inversion regularization analysis (Cracknell & Reading, 2015;Giraud et al., 2020;Tarabalka et al., 2009). However, to the best of our knowledge, the application of the topological relationships while deforming discrete units in geophysical inversion problems has not been addressed.
Among different orders of the topology (Thiele et al., 2016) we focus on first-order topology describing adjacency relationships between rock units. The utilized level-set inversion approach allows us to enforce small-scale topological rules as 185 morphological constraints to the inversion problem. To ensure that inverted models remain geologically realistic, we apply topological rules at each iteration. We take advantage of morphological rules of image processing techniques to prevent the nucleation of a given unit into another and for the model to obey topological rules. This becomes important for retaining the integrity of the predefined unit boundaries during the inversion and geological plausibility (age and deformation history).
Application of the mathematical morphology on geoscientific datasets have been evaluated based on the classification of input 190 data types (Heijmans, 1995;Serra, 1986;Soille and Pesaresi, 2002). However, it has rarely been applied within the geophysical inversion problems. Given that in a level-set inversion problem the model space can be assumed as a multinary image where the ordering of signed-distance values matters more than property values, morphological opening or closing (Vincent, 1993) rules can be applied on the multinary lithology ( ℎ ) at each iteration. For a 3D model, considering a structuring element of size ( × × ) 0< , , < , the closing notation (•) of ℎ can be shown as: 195 Where ⊝ and ⊕ demonstrate morphological dilation and erosion respectively (Jankowski, 2006;De Natale and Boato, 2017).
This operation ensures that all neighbouring model cells with the size of structuring element that does not fit in the background density, will be closed, thus nucleation of each lithology with size more than ( × × ) is prevented from occurring.

Synthetic examples 200
The aim of this section is to study the application of the introduced procedure on two synthetic case studies. We first benchmark the method using a well-known model and then simulate a more realistic application in a hard-rock scenario where 2D constraints are applied in a 3D inversion setting.

SEG/EAGE salt dome model
In this section, we use a simplified version of the SEG/EAGE salt dome model (Aminzadeh, 1996) in the same fashion as Li 205 et al (2016). This example demonstrates the application of the regularized level-set inversion for a model of two distinct rock units. We assume the salt dome in Fig Fig. 2c (from =3000, =0 to ′ =13000, ′ =11000 m) to image the salt body for the constraining purpose. We also add 5% noise with normal distribution to the forward calculated gravity datasets to generate the field datasets on the surface ( × = 30 × 30 = 900). It is assumed that the subsalt boundary and also the dipping part of the salt have been poorly imaged and are not interpretable from the seismic section (Fig. 2b). This restricts the constraining matrices from the seismic section only to the upper boundary of salt (Fig. 2d). The constraining matrices are then appended to the sensitivity 215 matrix as global and local regularizations for both units.
In many of the real case scenarios the starting model although follows the primary assumptions in the region, might be far from the reality. Thus, we start the inversion using a randomly generated starting model different from the true model to test and evaluate the similarity of the resulted geometry with the true geometry. In this example, we consider a disc as starting model ( Fig. 2e) that has intersections with two edges of the model boundaries. 220 The results of the level-set gravity inversion of the salt dome using starting model in Fig. 2e with and without the utilization of seismic information are illustrated in Fig. 3. The resulted model demonstrated in Fig. 3a is recovered after applying constraints along 2D seismic section based on the well-imaged upper boundary of the salt. We apply the constraints for this example along the seismic section using the weighting value (400) along the well-imaged top salt boundary (black cubes in

Example of hard-rock synthetic case
We have also tested the method on a second generated synthetic model (Fig. 4a) which contains four distinct geological bodies 245 with different density contrasts. The model simulates a hard-rock scenario and contains two exposed bodies (greenstones) with the same density contrast (330 . −3 ) surrounded by lower-density geological units (granitic background). This model is considered as the true model and the aim of the 3D generalized level-set inversion is to recover a model that is structurally close to the true model by constraining the inversion using the seismic information available in the area. The model is composed of × × = 40 × 30 × 40 = 48000 cubic model cells with 50 m resolution. A zero-offset straight synthetic seismic 250 section along y = 750 m using a finer grid mesh (10 m each cell dimension) is generated and 2% random noise with normal distribution is added to the amplitudes (Fig. 4c). The true model is then used to simulate gravity anomaly data at surface level ( × = 36 × 26 = 936) with 50 m spacing and assuming padding of 100 m. We add 5% noise with normal distribution to the gravity measurements (Fig. 4b).
As the main focus of this study is to test the application of 2D vertical constraints in a 3D level-set inversion, we only consider the seismic section for generating the starting model. For the purpose of testing our 255 algorithm, we simulate a starting model with an inaccurate representation of the geology of the area. In particular, connectivity of the greenstones and differences between the position and dipping of all layers (Fig. 4d). The starting model is used to initialize the signed-distance functions. To add distributed constraints along the seismic section to the inversion, we use the absolute value of reflectivity coefficients multiplied by constant values to directly translate boundaries from the seismic section to weighting matrices in the level-set formulation in Eq. (6). Translating these boundaries to weighting matrices has been 260 explained in detail in the theory section (Illustrated in Fig. 1b). Reflectivity as a measure of acoustic impedance can be calculated with the knowledge of the wavelet frequency. We use the same reflectivity matrix that was calculated for generating the synthetic seismic. Eventually, we apply weighting factors of 200 on model cells along the seismic section with sharp boundaries as maximum weights along reflectors for the constraining purpose. The morphological closing using a structuring element of size (100 × 100 × 100 3 ) is also applied to the resulting model to prevent nucleation of one lithology into the 265 other. The initialized value for is considered 35.05 . The inverted model stabilizes around an acceptable solution after 6 iterations with a total data misfit of 0.18 mGal (Fig. 5b). Qualitatively, visual inspection reveals that the resulting inverted model is in an acceptable agreement with the true model (Fig. 5a).
To assess the influence of seismic-derived constraints, we repeat the level-set inversion without applying the constraints along the section. The resulted model of 3D gravity inversion without seismic constraints is illustrated in Fig. 5d, with an overall 270 misfit of 0.10 mGal in data RMS error (Fig. 5e).
A comparison of final models from constrained and unconstrained inversion is presented in Fig. 6. The model differences between inverted models and the true model in Fig. 6c and Fig. 6d show more improvement in the constrained case. Model RMS error shows lower values for the constrained inversion (72.5 . −3 ). We also compare the results by measuring the structural similarity between the inverted models and the true model. A measurement for structural similarity (SSIM) (Wang 275 et al., 2004) that models the perceived changes in structural information of two different models A and B can be used as an indication of changes in the unit boundaries' location after level-set inversion.
( , ) = (2 + 1 )(2 + 2 ) � 2 + 2 + 1 �� 2 + 2 + 2 � (8) Where and are mean and variance of two models respectively, stands for the covariance of two models, 1 and 2 are two small constant values and is the dynamic range of the density contrast values. For two identical models, SSIM will become 1. The results indicate that final generated models from seismically constrained gravity inversion has recouped structural features of the true model up to 65%. Although the unconstrained inversion leads to a lower data misfit error, the structural similarity between the inverted model and true model is less favourable than in the seismically-constrained 285 case (33%). The difference between the structural similarities is an indication of the applicability of the approach to utilize

Geology of the area
The importance of Yilgarn Craton to the economy of Western Australia is evident as it consists of numerous granite-greenstone terranes hosting world-class deposits of gold and nickel (Whitaker, 2004). Several terranes are defined in this craton based on geochronological dating of magmatism and further geochemistry analysis (Cassidy et al., 2006;Pawley et al., 2007). The 305 eastern portion of the Yilgarn Craton is the Eastern Goldfields Superterrane which is divided into four terranes: Kalgoorlie, Kurnalpi, Burtville and Yamarna (Pawley et al., 2007). The Yamarna Terrane is located in the northeast (Fig. 7). It is proposed that the Yamarna Terrane evolved separately from the original Burtville terrane based on the different character of volcanosedimentary events (Pawley et al., 2012). While there are similarities between the Yamarna Terrane and the sequences around Kalgoorlie (a town host to one of the largest gold deposits in the world), the lack of significant historical mineral discovery is 310 at odds with the apparent prospectivity of the region. Complicated mineral exploration by extensive Phanerozoic cover (Lindsay et al., 2020), significant regolith (Anand and Paine, 2002) and poor outcrop encourage new geophysical investigations in this area to unlock the structure permissive for mineralizing processes (Goleby et al., 2004;Lindsay et al., 2020).  Numerous geophysical datasets including regional gravity data and 2D reflection seismic profiles have been collected in this area. However, few studies integrating these data have been done to estimate the 3D structure of the greenstones. The targeted greenstone belts are proximal to a series of major and minor faults and shear zones which make it necessary to utilize 320 geophysical integration for regional studies to have a better understanding of these metal hosts. Gravity datasets, although providing satisfying lateral resolution, are barely capable of estimating the depth of sub-horizontal interfaces. 2D seismic sections, on the other hand, contain information of deep structures that need complementary datasets to be laterally extended.

Geophysical datasets
The region of interest was chosen based on primary interpretation of greenstone locations approximated from the 2D seismic 325 section and surface geological mapping so that it covers an area where two major shear zones (Yamarna and Dorothy Hills) and bounded greenstones cross the long seismic section. This region of interest contains more sparse gravity measurements toward Burtville Terrane to the west and Dorothy Hills shear zone to the east (average 10km spacing) while denser gravity measurements with 2km spacing are collected around the Yamarna shear zone (Fig. 8). In this study, we use the interpolated regular grid of gravity anomaly data on the surface (from the Geological Survey of Western Australia) with 400-meter spacing 330

Starting model generation
A general overview of the geological setting of this region based on deep reflection seismic profiling with a focus on Laverton 340 Tectonic zone modelling is presented in Goleby et al (2004). We use their primary interpretation along the seismic section https://doi.org/10.5194/se-2021-65 Preprint. Discussion started: 16 June 2021 c Author(s) 2021. CC BY 4.0 License. within our region of interest together with modifications based on the integrated multi-scale study presented in Lindsay et al (2020) to generate a starting model for the gravity inversion in the Yamarna terrane.
Referring to primary seismic interpretations from Goleby et al (2004) and Lindsay et al (2020), major shear zones hardly extend deeper than 8km in depth. Focusing on the upper 8 km, we remove the effects of the regional gravitational trend and 345 recalculate the residual Bouguer anomaly (Gallardo and Thebaud, 2012). Therefore, we generate a model grid from the surface down to 10km depth and discretize to 500 m resolution with cubic cells.
Due to the insufficient number of studies in the region, we mostly rely on the results from a multi-scale study presented in Lindsay et al (2020) to generate the starting model. In this report, diverse geophysical datasets are utilized to conclude a general integrated interpretation of the region. The presented results from magnetotelluric modelling, reflection seismic reprocessed 350 image and potential field gravity inversion are used for geoscientific investigation. These results from the geophysical modelling are then used to relate petrophysical signatures to geological units using machine learning techniques. We initialize the density contrasts (with the base density of 2670 . −3 ) based on a primary gravity anomaly interpretation provided in the report as shown in Fig. 9a. The density contrast values assigned to defined units are inferred from the same report as well as interpretations from Goleby et al (2004) and summarized in Table 1 (first row). For more details about the integrated 355 interpretation of the region, we refer the readers to Lindsay et al (2020).
The ultimate aim of this case study is to produce a 3D model from the seismically constrained level-set inversion that is consistent with the surface gravity anomaly while including detectable, thus plausible, structures from the 2D seismic section.
To achieve this ambitious aim, we first break down the workflow to 3D property inversion followed by the unconstrained and constrained 3D level-set inversions. 360

Modifying the starting model
We use the generated knowledge-driven model as the prior model for the density contrast inversion. The Tomofast-x inversion platform (Giraud et al., 2021b;Giraud, et al., 2020;Giraud, et al., 2019a;Martin et al., 2020) is used for the property gravity inversion with smallness constraints from the prior model. The inverted density contrast model and the resulting calculated datasets together with the seismic overlayed image are shown in Fig. 9b to 9d. 365 Before level-set inversion, we modified the starting model based on the property inversion results and utilized the resulting model (Fig. 10) for the 3D level-set inversion. The main motivation for modifying the starting model based on the property inversion result is to generate a subsurface model that follows both geophysical datasets and also complements the prior geological knowledge of the region. We tie all this information together and obtain a geologically and geophysically plausible updated model of the area. In this modified model, we assign the same density contrast values to Mount Venn and Yamarna 370 greenstones and consider them as a single density contrast unit as suggested from the property inversion results. In addition, we propose an alternative scenario where we added another density contrast unit to the west as derived from the property https://doi.org/10.5194/se-2021-65 Preprint. Discussion started: 16 June 2021 c Author(s) 2021. CC BY 4.0 License. inversion result. Primary 3D level-set inversion of the region of interest using the initially defined starting model strongly suggests the generation of a new geological unit to the west with different dipping from the general east-dipping structure of the area. This appended unit to the west can also be inferred from the seismic image (Fig. 9d) where there are some flat-dipping 375 heterogeneities toward the west edge. We assign a data-driven density contrast heterogeneity to the unit toward the west of the model based on the property inversion results. The ultimate density contrast values used for generating the modified starting model are summarized in Table 1. This modified model is then used as the starting model for the 3D level-set inversion of the region forming five distinct level-set functions. The traversed seismic section in the region of interest covers all five defined rock units in the region (Fig. 10b). Furthermore, the modified starting model along this line is utilized for constraining purposes. 380 Out of convenience, we use the section view of the model along the seismic profile for most of the visualization purposes.
The discussion over testing different ranges of density contrast values for different rock units is beyond the scope of this study.
However, several sets of density contrasts based on generating a similar range of forward gravity datasets with the field datasets have been tested. Due to recovering almost the same structures by using different density contrast ranges, we present the results of one set of the resulting geometries (using densities in the second row of

Level-set inversion in 3D
The aim of this section is to illustrate the application of the gravity level-set inversion algorithm to constrain the 3D density model of the area with the seismic section. The starting model for the 3D level-set has five distinct units with densities as shown in Fig. 10. The model size for the 3D inversion is × × = 180 × 80 × 20 = 288,000, and the data size 395 is × = 50 × 20 = 1,000. As stated earlier in this section, very little is known about the area and the generated starting models are highly uncertain. Therefore, the results presented in this section are also undetermined and need to be interpreted with caution. The calculated gravity anomaly responses from the starting model and a comparison with the field datasets are shown in Fig. 11. A very high misfit between the observed and field datasets indicates that existing interpretations are far from reality.

Unconstrained 3D level-set inversion
We first implement the level-set inversion without applying constraints along the seismic profile to compare the resulting model from level-set inversions without and with seismic information. In this scenario, uniform local and global regularizations vectors are appended to the sensitivity matrix without including the weighting from seismic information. We set all 410 regularization parameters to 1 except for the layer at the surface level where we assign a small number (500) to the weighting terms to include information from geological mapping. Morphological closing is also applied as previously introduced using a structuring element of size (400 × 400 × 400 3 ) . The depth weighting term for compensating the decay of the gravitational force with depth is also added to the inversion problem as introduced and utilized in Boulanger & Chouteau (2001). The data RMS error of the unconstrained level-set inversion decreases smoothly converging at 1 mGal after 13 415 iterations (Fig. 13a). The small number of iterations for reaching the optimum value of the data RMS error is due to choosing a relatively broad interface between rock units ( = 2100 ). We previously discussed the importance of this parameter to control the search space of the problem in the theory section. Choosing such a large value for allows converging rapidly given the size of the model. The geometry of the density units of the resulted model however is not in agreement with the existing geological map and 420 seismic interpretations. The final modelled geometry of the greenstones from the unconstrained inversion is shown in Fig. 12a.
These geometries are terminated near the surface and are not extended to depth. This is mainly because we do not apply the constraints at depth. Therefore, the application of the regularizations from the seismic profile to improve the integrity of the units as they evolve laterally (as observed in the interpolated gravity grid) and vertically (as interpreted from the seismic profile) is necessary. 425

Constrained 3D level-set inversion
Next, we apply the spatial constraints along the seismic section to constrain the gravity inversion with the primary interpretation along seismic. We aim at reconciling gravity inversion with seismic interpretation. We use the same morphological constraints and depth weighting as we set for the unconstrained inversion in order to provide a fair comparison with the unconstrained inversion result. We also use the same value for the parameter. The updated regularization terms with weighted constraints 430 from the seismic information along the vertically extended profile are then applied to constrain the inversion. All arrays of vector are set as one except along the seismic section and for the layer at the surface level which we set as 2000 and 500 respectively, to suppress the changes of the signed-distance function and to constrain the inversion problem with spatially distributed regularizations. The weighting constants are tuned manually during the level-set inversion. As a rule of thumb, these constants are chosen between double to fourth of the cell size. Repeatedly, starting from a significant data RMS (more 435 than 20 mGal), it takes 12 iterations (Fig. 13b) for the data RMS error to reaches the minimum value (1.4 mGal). The final modelled geometry of the greenstones from the constrained inversion is also shown in Fig. 12b.

Figure 13: Evolution of data RMS error for the unconstrained (a) and constrained (b) 3D level-set inversion of Yamarna Terrane.
Comparing the final model from the inversion (Fig. 12b) and the prior model (Fig. 10), we can observe that the geometry of greenstones along the seismic profile is well-constrained and the surface features also align with the observed gravity datasets (Fig. 14). Most parts of the inverted model are in line with the existing assumptions in the area and represent an updated 445 geometry for greenstone units adjacent to the shear zones. There are differences between these two models along the seismic section while the general trend of features away from the seismic profile for both constrained and unconstrained scenarios show very similar patterns. This is an indication of the capability of the approach for applying the constraints along a vertical section and constraining the gravity inversion with sparsely distributed information within the model.   Fig. 15 the recovered geometry of the greenstones can be 455 correlated with some detectible features in the seismic image which were not primarily specified. We also provide a comparison of the recovered model with the integrated interpretations provided in Lindsay et al (2020) in Fig. 15b.
The source of the difference between these interpretations is mainly because in Goleby et al (2004) the focus was more on the whole of crust, whereas Lindsay et al (2020) were looking nearer to the surface and also utilized the re-processed seismic profile in their study.

Discussion
The main focus of this study has been to enable appending spatially distributed constraints in 2D from seismic section to 3D gravity inversion. We have shown this capability across two different synthetic models and a case study to demonstrate different workflows for constraining purposes depending on the available datasets and the area of study.
As detailed in the introduction, seismically constrained gravity inversion has been addressed for a long time. However, 470 constraining gravity inversion quantitatively with sparse seismic information in a cooperative workflow using the level-set technique is novel. While we have used information from reflection seismic as spatially distributed constraints, we see no limitation of using the same approach with other geophysical techniques that provide spatially distributed constraints within the model. For instance, any boundary recognition of rock units from passive seismic, magnetotelluric or electromagnetic techniques can also be translated to weighted constraints to the level-set inversion in the same fashion as we used for reflection 475 seismic. The level-set method is applied to a part of the eastern Yilgarn craton, a real-world case study where very little is known about the larger-scale 3D geology of the area. Testing revealed that the constraints generated from the 2D seismic profile and applied 490 to the 3D level-set inversion were effective and provided an informative update to the model of the subsurface. The results do not completely disagree with existing crustal interpretations presented in the literature (Goleby et al., 2004;Lindsay et al., 2020). We note that, based on our results, some new and modified features could be appended to complement existing interpretations.
From CDP number 13500 to 15000 to the left of the section in Fig. 15, the granite-gneiss outcrop is the main detectible feature 495 on the surface. However, we have assigned a small density contrast value (based on a small gravity anomaly) which also agrees with the moderate east-dipping structures in the seismic image (Unit A in Fig. 15b). The main detectible structure from the seismic image is almost compatible with the geometry of the rock units assigned to Yamarna Greenstone Belt (unit C) which in our resulting model is extended by unit D (a unit with the same density contrast as assigned to Dorothy Hills). This almost matches with the steep east-dipping structure defined in zone V from Lindsay et al (2020). Toward the right side of the section 500 after Yamarna Greenstone Belt (from CDP: 16700 to 18000 and from 5 to 10 km), there have been minor interpretations that could explain the flat-lying structures of the seismic image. In the resulting model, these structures are assigned to unit F and unit G with the same density contrast assigned to Dorothy Hills and Lake Yeo respectively.
The main detectible feature from the constrained level-set inversion as presented in Fig. 12b is the extension of the density contrast that we assign to the Yamarna Greenstone Belt unit that displays a keel-like shape with a dip to the southeast. The 505 western edge of the feature also results in strong reflectors in the seismic image (Units C and D in Fig. 15) and has been interpreted to be the Yamarna Shear Zone (Goleby et al., 2004;Lindsay et al., 2020). Our result indicates that the predicted and modelled Yamarna shear zone aligns with density contrast patterns that follow anomalously high gravity response laterally https://doi.org/10.5194/se-2021-65 Preprint. Discussion started: 16 June 2021 c Author(s) 2021. CC BY 4.0 License. away from the seismic survey trace. The results suggest that the greenstone belt extends, along with the shear zone, up to 8 km to the north and south and likely beyond the geographic boundaries of the model. 510 The other interesting feature from the inverted model also in strong agreement with surface geological maps and geophysical information is the depth extension of the Dorothy Hills Greenstone Belt (Units E). Due to the narrow width of this greenstone belt on the surface, the extension of this unit toward deep parts of the model is unlikely. Referring to studies that address greenstone structures in Canada (Thurston, 2015) and Australia (Blewett et al., 2010;Gallardo & Thebaud, 2012) the assumption is that the width of the greenstone belts are indicative of their depth which explains the shallow depth of unit E. 515 A disconnected extension of this unit toward the eastern part of the model together with the density contrast that we assigned to the Lake Yeo unit have formed volumes with high density contrast (Units F and G). This also supports the existing interpretation provided in Lindsay et al (2020) where they assign a higher density domain to the eastern side of the Dorothy Hills shear zone.
While our level-set inversion results strongly suggest the creation of a density contrast unit with a different dipping structure 520 (Unit A) to the west of the model there is a lack of evidence in the literature about the extent of such a density contrast unit. This is mainly because toward the west edge of the model there is a granite-gneiss outcrop and surface geological evidence ( Fig. 7, toward Burtville Terrane) fails to confirm the existence of such density contrast in depth. Knowledge about physical properties in a hard-rock environment is not enough to constrain the geological concepts and the effects of other processes should also be considered (Dentith et al., 2020). This recovered unit, be it a data-driven unit from observed gravity datasets 525 lies in an area with lower density contrast and higher magnetic susceptibility domain based on presented results of gravity and magnetic inversion in Lindsay et al., (2020). These changes in density and magnetic susceptibility might not necessarily support the introduction of new geological unit and could be an indication of some secondary geological processes (Dentith et al., 2020;Saltus and Blakely, 2011;Whitaker, 2004) that results in local bulk heterogeneities regarding the complex mineral compositions around this region. 530 Opportunities for improving the presented level-set procedure are centred on interpretation and model uncertainty. Firstly, there is considerable uncertainty regarding the appropriate velocity model used during the seismic imaging process utilized within the level-set algorithm. Some existing potential field geophysical inversion methods (Giraud et al., 2019) use uncertainty as a constraint, and integration with reflection seismic data would help resolve more complicated scenarios involving sparse data, limited petrophysical contrast and inadequate model assumptions. In structurally complex areas such as 535 in hard-rock scenarios, there are high uncertainties regarding the local and regional interpretations as well as in seismic imaging techniques. This indicates the necessity of quantified inclusion of uncertainty regarding the interpretation and imaging process of seismic data while constraining the level-set inversion of potential field datasets. It suggests the development of a new technique that enables interactive constraining of the potential field and reflection seismic datasets within a level-set algorithm.
One advantage of enabling this type of integration in the level-set framework is the compatibility of the results with implicit geological modelling. Seismic interpretations and gravity datasets are both very common datasets for generating 3D models of an area while there is no technique that allows including seismic uncertainty within 3D integrated implicit modelling algorithms.
The presented results on different datasets revealed that constraining level-set inversion with sparse constraints from lowuncertainty datasets can effectively improve the (geological) plausibility of the results. Adding constraints to the level-set 545 inversion seemed to be effective in all scenarios and is even more obvious with increasing complexities in the models. Sparse constraints, which we see no limitation for their source, force the inverted models to retain the available prior knowledge of the area to a high degree. This introduces a resulting model incorporated from several sources of geological and geophysical information and is a step toward automating the process of 3D geo-modelling.

Conclusion 550
We have presented the utilization and extension of a generalized level-set approach for seismically-constrained gravity inversion across different scenarios. The flexibility of the level-set approach we followed allowed us to append 2D constraints from spatially distributed seismic information to 3D level-set gravity inversions. We tested the method on two different synthetic models with different scenarios and levels of complexity. The proposed technique was then applied to the geologically complex Yamarna terrane and an updated model of the area with supporting discussions was provided. The 555 proposed approach has proven to be reliable to quantitatively include sparse constraints with low uncertainty to the level-set inversion. In addition, there is considerable flexibility regarding the constraints which makes the approach widely applicable for integrating with other geophysical datasets and geological models. The availability of other sources of information in the studied area about the depth of interfaces in the future can also be utilized in this framework for possible modifications of the recovered model. Also, any improvement in the interpretation of the seismic profile based on future evidence can be directly 560 injected into the inversion to improve the final recovered model.

Competing interests
The authors declare that they have no conflict of interest. MR designed the methodology and performed all modelling and inversion presented in the manuscript with the continuous assistant of JG. MR is the main writer of the manuscript which was redacted with the support of all co-authors. ML assisted in 570 the generation of the 3d geological models of Yamarna terrane and interpretation of the results. MJ provided guidance and supervision while the project was being carried out. VO assisted with the mathematical notations in the methodology section.
The manuscript is part of the MR PhD project supervised by the rest of the authors.

Acknowledgment
The work has been supported by the Mineral Exploration Cooperative Research Centre whose activities are funded by the 575 Australian Government's Cooperative Research Centre Program. This is MinEx CRC Document 2021/xx.