Utilisation of probabilistic MT inversions to constrain magnetic data inversion: proof-of-concept and field application

. We propose, test and apply a methodology integrating 1D magnetotelluric (MT) and magnetic data inversion, with a focus on the characterization of the cover-basement interface. It consists of a cooperative inversion workflow relying on standalone inversion codes. Probabilistic information about the presence of rock 20 units is derived from MT and passed on to magnetic inversion through constraints combining such structural constraints with petrophysical prior information. First, we perform the 1D probabilistic inversion of MT data for all sites and recover the respective probabilities of observing the cover-basement interface, which we interpolate to the rest of the study area. We then calculate the probabilities of observing the different rock units and partition the model into domains defined by combinations of rock units with non-zero probabilities. Third, we combine 25 such domains with petrophysical information to apply spatially-varying, disjoint interval bound constraints to least-squares magnetic data inversion. We demonstrate the proof-of-concept using a realistic synthetic model reproducing features from the Mansfield area (Victoria, Australia) using a series of uncertainty indicators. We then apply the workflow to field data from the prospective mining region of Cloncurry (Queensland, Australia). Results indicate that our integration methodology efficiently leverages the complementarity between separate MT 30 and magnetic data modelling approaches and can improve our capability to image the cover-basement interface. In the field application case, our findings also suggest that the proposed workflow may be useful to refine existing geological interpretations and to infer lateral variations within the basement. expected, which are formulated on a case 150 by case basis. We first calculate the depth to basement interface probability 𝑝 𝑖𝑛𝑡 for each MT sounding. These probabilities are derived for each MT sites and can be interpolated on the mesh used for magnetic inversion. The interpolation of such probability distributions can be performed using different approaches and integrate various types of geophysical or geological constraints. In this work, we use a linear interpolation scheme in the synthetic case study, and the Bayesian Estimate Fusion algorithm of Visser and Markov (2019) in the case study. We note 155 that other techniques could be used for a similar purpose, such as the Bayesian Ensemble Fusion (Visser, 2019; Visser et al., 2021), or discrete and polynomial trend interpolations (Grose et al., 2021). Following this, assuming a sedimentary basin lying on top of the basement, we can define for each model cell the probability of being in presence of the basement 𝑝 𝑏𝑠𝑚𝑡 as: et al., 2008; Lajaunie et al., 1997). It is constituted of a sedimentary syncline abutting a faulted contact with a folded basement. The model’s complexity was increased with the addition of a fault and an ultramafic intrusion. Details about the original 3D geological model are provided in Pakyuz-Charrier 235 et al. (2018b). Here, we increased the maximum depth of the model to 3150 m and added padding in both horizontal directions. Figure shows the non-padded 2D section extracted from the reference 3D geological model. a similar workflow as presented in sections 2.1 and 3.3, they derived a cover-basement interface distribution for each MT site using the model posterior distribution obtained using 1D probabilistic inversions. Then, a 395


Introduction
Geophysical integration has been gaining traction in recent years, be it from the joint or cooperative point of 35 views. A number of approaches for the joint modelling have been developed with the goal of exploiting the complementarities between different datasets (see for instance the reviews of Farquharson, 2016, andMoorkamp et al., 2016, andreferences therein). As summarized in the review of Ren and Kalscheuer (2019), "joint inversion of multiple geophysical datasets can significantly reduce uncertainty and improve resolution of the resulting models", be it for the modelling of a single property (e.g., resistivity for joint controlled-source 40 electromagnetic and magnetotelluric MT data, or density for joint gravity anomaly and gradiometry data), or of multiple properties (e.g., the joint inversion of seismic and gravity data to model P-velocity and density). In the second case, joint inversion approaches can be grouped into two main categories based on the hypothesis they rely on. Structural approaches allow to jointly invert datasets with differing sensitivities to the properties of the subsurface through the premise that geology is such that spatial changes in inverted properties should be 45 collocated. Structural constraints can then be used as a way to link two or more datasets jointly inverted for by encouraging structural similarity between the inverted models (Haber and Oldenburg, 1997;Gallardo and Meju, 2003). Alternatively, petrophysical approaches utilise prior petrophysical information (e.g., from outcrops, boreholes, or the literature) to enforce certain statistics in the recovered model so that it resembles the petrophysical measurements' (Lelièvre et al., 2012;Sun and Li, 2015;Giraud et al., 2017;Astic and Oldenburg, 50 2019). Whereas structural and petrophysical approaches are well suited to exploit complementarities between datasets in a quantitative manner, running joint inversion might be, in practice, challenging and requires significantly more computing power than the separate inversions.
In this contribution, we present a new multidisciplinary modelling workflow that, instead of the more common simultaneous joint inversions, relies on sequential, cooperative modelling. It follows the same objectives as the 55 two categories of joint inversion mentioned above in that structural information is passed from one domain to the other and it uses petrophysical information to link domains. The development of the sequential inversion scheme we present was motivated by a similar idea as Lines et al. (1988) who states that "the inversion for a particular data set provides the input or initial model estimate for the inversion of a second data set". A further motivation was to design a workflow capable of integrating the inversion of two or more datasets quantitatively using 60 standalone modelling engines that run independently. The general concept behind the methodology we present is summarised in Figure 1.
In this paper, the workflow is applied to the sequential inversion of magnetotelluric (MT) followed by magnetic data, taking into account the importance of robustly constraining the thickness of the regolith in hard rock imaging and mineral exploration. This is motivated by the relative paucity of works considering cooperative workflows to 65 integrate MT and magnetic data together with the recent surge in interest for the characterisation of the depth-tobasement interface in mineral exploration, despite these two geophysical methods being part of the geoscientists' toolkit for depth-to-basement imaging. Historically, MT has often been integrated with other electromagnetic methods or with seismic data, and with gravity to a lower extent (see review of the topic of Moorkamp, 2017). It is, however, seldom modelled jointly with magnetic data unless a third dataset is considered (e.g., Oliver-Ocaño 70 et al., 2019 ;Zhang et al., 2020;Gallardo et al., 2012;Le Pape et al., 2017). We surmise that this is because (1) the interest for integrating MT with other disciplines arose primarily in oil and gas and geothermal studies and https://doi.org/10.5194/se-2021-124 Preprint. Discussion started: 3 November 2021 c Author(s) 2021. CC BY 4.0 License. relied on structural similarity constraints for reservoir or (sub)salt imaging, (2) of the difference in terms of spatial coverage between the two methods elsewhere, (3) the differences in terms of sensitivity to exploration targets and (4) the difficulty to robustly correlate electrical conductivity and magnetic susceptibility. Bearing these 75 considerations in mind, we developed a workflow incorporating MT and magnetic inversion with petrophysical information and geological prior knowledge.
In the workflow we develop, we exploit the differences in sensitivity between MT and magnetic data. On the one hand, the MT method, used in a 1D probabilistic workflow as presented here, is well-suited to recover vertical resistivity variations and interfaces, especially in a sedimentary basin environment . MT 80 data are, however, poorly sensitive to resistors, especially when they are overlaid by conductors (e.g., , which makes it difficult to differentiate between highly resistive features, such as intra-basement resistive intrusions. On the other hand, magnetic data inversion is more sensitive to lateral magnetic susceptibility changes and to the presence of vertical or tilted structures or anomalies. Bearing this in mind, we first derive structural information across the studied area in the form of probability distributions of the interfaces between 85 geological units, extracted from the interpolation of the results the 1D probabilistic inversion of MT data. From there, the probability of occurrence of geological units can be estimated in 2D or 3D. These probabilities are used to divide the area into domains where only specific units can be observed (e.g., basement, sedimentary cover, or both). Such domains are then passed to magnetic data inversion, where they are combined with prior petrophysical information to derive spatially varying bound constraints that are enforced using the alternating direction method 90 of multipliers or 2D or 3D inversion (see Ogarko et al., 2021a, for application to gravity data using geological prior information). Finally, uncertainty analysis of the recovered magnetic susceptibility model is performed and rock unit differentiation allows to control the compatibility of magnetic inversion results with the MT data. The workflow is summarised in Figure 1. The remained of this paper is organised as follows. We first introduce the methodology and summarise the MT and magnetic standalone modelling procedures we rely on. We then introduce the proof-of-concept in detail using a realistic synthetic case study based on a model of the Mansfield area (Victoria, Australia), which we use to explore the different possibilities for integrating MT-derived information and petrophysics offered by our 100 workflow. Following this, we present a field application using data from Cloncurry (Queensland, Australia) where we tune our approach to the specificity of the area. Finally, this work is placed in the broader context of geoscientific modelling and perspectives for future work are exposed in the discussion section.

MT inversion for interface probability 105
The MT method is a natural source electromagnetic method. Simultaneous measurements of the fluctuations of the magnetic and electric fields are recorded at the Earth's surface under the assumption of a source plane wave.
The relationship between the input magnetic field H and the induced electrical field E, which depends on the distribution of the electrical conductivities in the subsurface, is described by the impedance tensor Z, as follows: Resistivity models derived from MT data are found by forward modelling and inversion of the impedance tensor 110 Z, generally using gradient-based deterministic methods (see, e.g., Rodi et al., 2012). Deterministic approaches provide a single solution which minimizes the objective function considered during the inversion, but limited information of the uncertainty around this model can be derived. A global characterization of the uncertainty is possible using a Bayesian inversion framework, but its expensive computing cost limits its application to approximated and/or fast forward modelling solvers (Conway et al., 2018;Manassero et al., 2020;Scalzo et al., 115 2019). In this study we alleviate this by considering a 1D behaviour of the Earth. While this assumption can stand within layered sedimentary basins, it may fail in more complex geological environments. As we describe below, we account for this source of uncertainty within our Bayesian inversions.
Within the context of Bayesian inversion, the solution of the inverse problem consists in a posterior probability distribution, calculated from an ensemble of models that fits the data within its error. The posterior probability 120 distribution p( | ) is obtained using Bayes' theorem, defined as: The prior distribution p( ) contains prior information on the model parameters . In this work, we assume a relatively uninformed prior knowledge, using a uniform prior distribution on the electrical resistivity bounds with values set between 10 −2 and 10 6 Ω. m. Using a uniform prior with such wide boundaries allows the inversion to be mainly data-driven and to remain independent from assumptions about the distribution of electrical 125 resistivity into the Earth. A Gaussian likelihood function p( | ) defining the data fit is used is used: The term inside the exponential is the data misfit, which is the distance between observed data and simulated data ( ), scaled by the data covariance matrix , which defines data errors and their correlation across frequencies. We consider two main sources of uncertainty to calculate :  data processing errors, which we model introducing a matrix ; 130  errors introduced by the violation of the 1D assumption when using 1D models, which we model introducing .
To calculate , we first define , the covariance matrix accounting for the EM noise and measurement errors, which is estimated during MT data processing. In this study we assume uncorrelated processing noise across frequencies, reducing to a diagonal matrix. Following this, we define as the covariance matrix accounting 135 for the discrepancy between 1D models and the multi-dimensional Earth the data is sensitive to. The analysis of the characteristics of the MT phase tensor allows to detect 2D and 3D effects present in the MT data as a function of frequency (see Caldwell et. al., 2004).
Following a dimensionality error model developed by  for depth to cover mapping applications, we translate the phase tensor parameters into dimensionality uncertainties to compensate for the 140 limitations of the 1D assumption when performing 1D inversion. This dimensionality uncertainty is added to the existing data uncertainty , such that the inversion considers both sources of uncertainty in , which we calculate as = + . The 1D MT trans-dimensional Markov chain Monte Carlo algorithm used in this study  solves for the resistivity distribution at depth, and the number of layers in the model and takes as input the determinant of the impedance tensor Z. The output of the probabilistic inversion 145 consists in an ensemble of models describing the posterior probability distribution.
In this paper, we are interested in the depth of the basement. For each MT sounding, each model of the model ensemble is analysed and the transitions from a conductive sedimentary layer into a resistive basement is extracted in order to construct a probability distribution of the depth to basement interface. This feature extraction relies on assumptions made on the electrical resistivity of the different lithologies expected, which are formulated on a case 150 by case basis. We first calculate the depth to basement interface probability for each MT sounding. These probabilities are derived for each MT sites and can be interpolated on the mesh used for magnetic inversion. The interpolation of such probability distributions can be performed using different approaches and integrate various types of geophysical or geological constraints. In this work, we use a linear interpolation scheme in the synthetic case study, and the Bayesian Estimate Fusion algorithm of Visser and Markov (2019) in the case study. We note 155 that other techniques could be used for a similar purpose, such as the Bayesian Ensemble Fusion (Visser, 2019;Visser et al., 2021), or discrete and polynomial trend interpolations (Grose et al., 2021). Following this, assuming a sedimentary basin lying on top of the basement, we can define for each model cell the probability of being in presence of the basement as: with the cumulative distribution function of the probability function , from the surface downwards. 160 Consequently, for each cell within the model, the probability of being in presence of sedimentary rocks, , is given as = 1 − . This allows us to derive three domains characterised by = 0 (basement only), ∈ ]0,1[ (basement and non-basement units allowed) and = 1 (non-basement only) that will define the intervals used for the bounds constraints applied to magnetic data inversion as summarised below. 165

Formulation of the magnetic data inverse problem
In this section, we summarise the method used to enforce disjoint interval bound constraints during magnetic data inversion. We largely follow Ogarko et al. (2021a), which we extend to locally weighted bound constraints. The geophysical inverse problem is formulated in the least-squares sense (see chap. 3 in Tarantola, 2005). The cost function we minimize during inversion is given as: 170 where is the observed data and ( ) the forward response produced by model , a vector of ℝ , with is the total number of model-cells. The second term corresponds to the model damping (or smallness) term, with weight α , depth weighting operator ; is the prior model. The operator is defined following Portniaguine and Zhdanov (2002). The third term is the gradient damping (or smoothness) term.
The diagonal matrix adjusts the strength of the regularization. It is modulated by α (introduced in details in 175 the next section).
We solve eq. (4) while constraining the inversion using multiple bound constraints. The problem can be expressed in its generic form as: where ℬ is the interval, or set of intervals, binding the i th model-cell. The general definition of is: , with , > , , ∀ ∈ [1, ] and ∈ ℐ where , and , define, the bounds of the inverted property is the index of the rock unit; is the number of 180 bounds in the i th interval; ℐ contains the indices of the model-cells where the bound constraints are applied. A summary of the algorithm solving this problem using the alternating direction method of multipliers technique (ADMM) is given in Appendix A.

Integration with MT modelling
The set of intervals ℬ from eq. (6) and (7) can be defined homogenously across the entire model or accordingly 185 with prior information. In the latter case, it allows to define spatially-varying bound constraints and to activate them only is selected parts of the study area. In the case presented by Ogarko et al. (2021a), probabilistic geological modelling was used to determine such bounds constraints for gravity inversion. The approach we propose here follows the same philosophy. Instead, we use probabilistic MT modelling, which can be used to where , is the observation probability for the th rock unit; , is a threshold value above which the probability is considered sufficiently high to be considered. In what follows we use , = 0.

Model misfit 195
In the synthetic study we present, we evaluate the capability of inversions to recover the causative model using the root-mean-square (RMS) misfit between the true and inverted models (RMS model misfit). We calculate this indicator as: where and are, respectively, the true and inverted models.

Membership analysis 200
In the context geophysical inverse modelling, membership analyses provide a quantitative estimation of interpretation uncertainty to interpretation of recovered petrophysical properties. Using the intervals defined in eq. (7) and eq. (8), we apply a triangular function to the inverted magnetic susceptibilities to determine the membership values to the different rock types. For the i th model cell, the membership value to the j th rock type ω is calculated using ramp functions as: 205

Total model entropy
Using the membership values ω, we calculate the total model entropy of the model, , which is the arithmetic mean of the information entropy (Shannon, 1948) of all model-cells: which is as a measure of geological uncertainty in probabilistic models and of the fuzziness of the interfaces (Wellmann and Regenauer-Lieb, 2012). 210 https://doi.org/10.5194/se-2021-124 Preprint. Discussion started: 3 November 2021 c Author(s) 2021. CC BY 4.0 License.

Jaccard distance
In addition to calculating , the membership values ω can be used to interpret the inversion results in terms of rock units. The index k of the rock unit the model-cell with a given inverted magnetic susceptibility value can be interpreted as is given as, for the i th model-cells: Calculating the index of the corresponding rock unit in each model-cell, we obtain a rock unit model . 215 Using and (the latter being the true rock unit model), we calculate the Jaccard distance (Jaccard, 1901), which is a metric quantifying the similarity between discrete models. In the context of geological modelling, it is reflective of the dissimilarity between geological models (Schweizer et al., 2017). Here, we use it to compare the inverted model and the true model as follows: where ⋂ and ⋃ are the intersection and union of sets, respectively; |•| is the cardinality operator, measuring the 220 number of elements satisfying the condition. A useful and simple interpretation of J is that it represents the fraction of cells assigned with the incorrect rock unit.

Synthetic case study
The synthetic case study that we use to test our workflow was built using a structural geological framework initially introduced in Pakyuz-Charrier et al. (2018. It presents realistic geological features that reproduce field 225 geological measurements from the Mansfield area (Victoria, Australia). The choice of resistivity and magnetic susceptibility values to populate the structural model was made to test the limits of this sequential, cooperative workflow and to show its potential to alleviate some of the limitations inherent to potential field and MT inversions. To this end, we selected a part of the synthetic model which features violate the hypothesis of a 2D model to challenge the workflow we propose. 230

Survey setup
The structural geological model was derived from foliations and contact points using the Geomodeller® software Guillen et al., 2008;Lajaunie et al., 1997).  Figure 2a shows the non-padded 2D section extracted from the reference 3D geological model.
We assigned magnetic susceptibility in the model considering non-magnetic sedimentary rocks in the basin units (lithologies 3, 5 and 6 in Table 1) and literature values (see Lampinen et al., 2016) to dolerite (lithology 4), diorite 240 (lithology 2) and ultramafic rocks (lithology 1). We assign electrical resistivities assuming relatively conductive https://doi.org/10.5194/se-2021-124 Preprint. Discussion started: 3 November 2021 c Author(s) 2021. CC BY 4.0 License. sedimentary rocks and resistive basement and intrusive formations. Resistivities in sedimentary rocks might vary orders of magnitude, and mainly depend on porosity, which is linked to the degree of compaction and the type of lithology, and the salinity of pore fluid (Evans et al., 2012). The three sedimentary layers were assigned different resistivities values, of 30 Ωm, 10 Ωm and 50 Ωm for basin fill 3, 2 and 1 respectively (see Table 1), basin 1 being 245 the oldest and deeper formation. Metamorphic and intrusive rocks as found in the crust generally present high resistivities (Evans et al., 2012). In what follows we focus on the modelling of data located along the line shown in Figure 2, simulating the modelling magnetic data in 2D, while considering 3D MT data. The modelled rock units and their petrophysical properties are given in Table 1. The geological, magnetic susceptibility and resistivity true models are shown in Figure 2. 250    Figure 2.
Airborne magnetic data were simulated for a fixed wing aircraft flying at an altitude of 100 m above topography.
The forward solver we use follows Bhattacharyya (1964) to model the total magnetic field anomaly. In this example, we model a magnetic field strength equal to 57,950 nT. It corresponds to the International Geomagnetic Reference Field for the Rawlinna station, Western Australia. 270 We add normally distributed noise with an amplitude equal to 2.5% of the average amplitude of the data. We simulate noise contamination by adding noise sampled randomly from by a normal distribution characterised by a standard deviation of 3.8 nT and a mean value of 0 nT. For the simulation of geological "noise", we then apply a Gaussian filter to such random noise to obtain spatially correlated values. The uncontaminated and noisy data are shown in Figure 3. For the inversion, the objective data misfit is set accordingly with the estimated noise. 275

MT data
The synthetic MT data is computed using the complete 3D resistivity model derived from the 3D geological 280 model, to simulate the influence of 3D effects. The core of the electrical conductivity model used the same discretization as the magnetic susceptibility model. Relationships between geological units and electrical resistivities follow Table 1. A padding zone was added in both horizontal directions and vertically, in order to satisfy the boundary conditions required by the forward solver. The final mesh has 160 x 160 x 62 cells in the , and directions, respectively. The ModEM 3D forward modelling code (Egbert and Kelbert, 2012;Kelbert et 285 al., 2014) Figure 2). The frequency range spans from 10 KHz to 0.01 Hz, with 6 frequencies per decades, for a total of 37 frequencies; 5% magnitude Gaussian white noise was added to the synthetic data before running the 1D inversions.
In the following subsections, we present the results of the modelling of synthetic data along a 2D section (see 290 Figure 2c) from the 3D resistivity volume, following the workflow proposed in Sect. 2. Along this section, 16 MT sites are used. We start with the modelling of MT data to derive constraints and prior information for the inversion of magnetic data.

1D Probabilistic inversion of MT data and derivation of cover-basement interface probabilities
We perform the 1D MT inversions of each MT sounding independently, using a 1D trans-dimensional Markov 295 chain Monte Carlo sampler . The presence of 2D and 3D effects was compensated for using larger uncertainties, calculated using the workflow introduced in Seillé and Visser (2020) and summarized in Sect. 2.1.
All the 1D trans-dimensional Bayesian inversions ran using 60 Markov chains with 10 6 iterations each. For each chain, a burn-in period of 750,000 samples (75% of the total) was applied to ensure convergence, after which we 300 recorded 100 models equidistantly spaced. The model ensembles are then constituted by 6000 models for each MT site. An example of the posterior distribution for three MT sites is shown in Figure 4. The distribution of interfaces within the posterior ensemble of 1D models is described by a change-point histogram. Figure 4a and From the posterior ensembles of models and interfaces, a depth to basement probability distribution is calculated 305 independently for each MT site. For this synthetic case, we simply assume that the transition from the sedimentary cover to the basement occurs when a layer L1 of resistivity 1 < is followed by a layer L2 of resistivity 2 > , with = 200 Ωm. As mentioned in Sect. 2.1, this process is applied to each model of the ensemble and allows the extractions of geological features of interest from the posterior model ensemble. If less than 0.1 % of the all transitions observed in the ensemble presents the feature defined earlier using , we then assume that the 310 transition is not observed. This situation occurs for MT sites MT14, MT15 and MT16 (see Figure 4c for site MT16). Figure 4 shows the interface probability and the depth to basement interface probability for three selected sites along the profile.

Bound constraints 325
Starting from values calculated for each MT site, we interpolate the probabilities from MT onto the mesh used for magnetic data inversion. In this synthetic example, we use a linear interpolation scheme. The interpolated probabilities from MT are shown in Figure 5.
The interpolated probabilities are used to define domains for the application of bound constraints during magnetic data inversion. In this example, we complement information from MT inversions with the assumption that dolerite 330 outcrops and is mapped accurately (unit 4, intrusive, see Table 1 and Figure 2). Using this, we adjust the domains for the corresponding few model-cells at surface level, only at two locations where we assume these outcrop to be known. The domains for the bound constraints we obtain are shown in Figure 5c, domains 1 and 2 indicate parts of the model where MT inversion suggests a single rock unit. This means that in the corresponding model cells, a single interval will be used in the definition of the bound constraints, while two intervals will be used 335 otherwise.

Prior model
The prior model for magnetic inversion is obtained using the MT derived rock unit probabilities (Figure 5a-c) and the magnetic susceptibility of the rock units given in Table 1. We calculate it as follows, for the ith model-cell: We chose to use , , the lower bound for each rock unit (or group of rock units), as it constitutes the most 340 conservative assumption about magnetic susceptibility from the range of plausible magnetic susceptibilities. The resulting prior model is shown in Figure 5d.

Inversion of magnetic data and uncertainty analysis
In this section, we study the influence of MT-derived prior information onto magnetic inversion and estimate the related reduction of interpretation uncertainty. The different scenarios tested here are summarized in Table 2, while the corresponding in version results are shown in Figure 6. 350   proposed workflow has the capability to improve the recovery of the sedimentary cover thickness significantly when compared to cases that do not use MT-derived ADMM constraints across the entire model (Figure 7a, b and   c). 370 To complete our analysis, we propose a simple metric combining two statistical indicators by weighting the RMS model misfit resulting from geophysical inversion ( ) with the corresponding entropy values ( ), which reflect geological interpretation uncertainty. Figure 8 shows the apparent relationship between the amount of prior information and constraints used in inversion and the proposed indicator. We observe a continuous reduction of × from case (a) to (f), indicating a progression towards models that are closer to the true model and less 375 ambiguous to interpret. Overall, together with values of the Jaccard distance J (Figure 6), our testing reveals that, using MT-derived ADMM constraints provides inversions results that are the closest to the true model. However, it also suggests that from case (c) to (f), the impact of additional constraints decreases comparatively to cases (a) to (c). It also indicates that using the prior model has a limited effect on inversion, as cases (a)

Field case study
We propose an application example illustrating potential utilisations of the proposed sequential inversion workflow. Using observations made in the synthetic case, our aim here is to integrate MT with magnetic data 385 inversion using the case relying on MT-derived ADMM bound constraints with homogenous starting model and smoothing constraints. We note that, to the best of our knowledge, this is the first time that the ADMM algorithm is used to constrain magnetic data inversion.
We used existing results of the depth to basement derived using MT within a probabilistic workflow (Seillé et al., 2021) in the Cloncurry area in Queensland (see Figure 9). These results are used to constrain the magnetic 390 inversion.

Geoscientific context and area of interest
The depth to basement interface probability used for the magnetic inversion was derived by Seillé et al. (2021).
Using a similar workflow as presented in sections 2.1 and 3.3, they derived a cover-basement interface distribution for each MT site using the model posterior distribution obtained using 1D probabilistic inversions. Then, a 395 https://doi.org/10.5194/se-2021-124 Preprint. Discussion started: 3 November 2021 c Author(s) 2021. CC BY 4.0 License. probabilistic depth to basement interface across the survey area was derived combining these MT depth to basement estimates with drill-hole depth-to-basement estimates, constrained by a structural model in the form of interpreted faults across the area. The combination of estimates coming from different sources of information to derive a posterior distribution on the depth to basement was performed using the Bayesian Estimate Fusion algorithm of Visser and Markov (2019). This probabilistic workflow allowed to obtain a probabilistic map of the 400 depth-to-basement interface. The structural model used during the Estimate Fusion assumed vertical faults, which is a valid assumption given the near vertical behaviour of the main faults in the area (Austin and Blenkinsop, 2008;Case et al., 2018).
The geological map of the area is shown in Figure 9. In this study, we focus on a 2D profile (L26, see location on map in Figure 9a), and invert the corresponding 410 magnetic data extracted from the anomaly map shown in Figure 10a and Figure 10b. The choice of an East-West oriented profile was motivated by the North-South orientation of the main structures in the area and by the geological features the known geology and the geophysical measurements suggest. The profile is nearly perpendicular to these structures, making it suitable for use within a 2D inversion scheme. It crosses the North-South oriented Mount Margaret Fault, which is thought to belong to the northern part of the regional Cloncurry 415 Fault structure, a major crustal boundary that runs North-South over the Mount Isa Province (Austin and Blenkinsop, 2008;Blenkinsop, 2008). This boundary separates two major Paleoproterozoic sedimentary sequences (Austin and Blenkinsop, 2008). The geological modelling performed by Dhnaram and Greenwood (2013)  the West and the Soldier Caps Domain to the East. In our study area, the Constantine domain is covered by non-420 magnetic cover constituted by Mesozoic and Cenozoic sediments, lying on what is believed to be consitituted by the Mount Fort Constantive volcanics, in some places intruded by the Williams supersuite pluton. On the eastern side, the Soldier Caps Domain is also covered by Mesozoic and Cenozoic sediments, and the basement is interpreted to be a succesion of volcanic and metamorphic rocks (Dhnaram and Greenwood, 2013). The depth to basement probabilistic surface derived by Seillé et al. (2021) along the W-E profile (see Figure 10c) 425 presents a relatively shallow basement in the western part of the profile (~ 100 meters thick, with some lateral variations). In the eastern part of the profile, the model indicates that a two-steps fault system controls the thickening of the basin to the east. It reaches ~ 350 meters thickness in the eastern part. The depth to basement model along the profile shown in Figure 10c is relatively well constrained by MT and drill hole data. However, the interpolation method we used imposed spatial continuity between estimates that, using MT soundings with 430 relatively large separation (2 km) and sparse drill holes, did not allow for the definition of small-scale depth to basement lateral variations. Significant lateral variations were allowed during the interpolation using the fault traces indicated by the structural data, defining areas where discontinuities are expected (in the location of the faults), to encourage a relaxation of the spatial continuity and to allow for sharp jumps. The magnetic data shown in Figure 10a suggest that other faults and other lateral discontinuities could exist. 435 In this work, we assume a non-magnetic sedimentary cover, and a magnetic basement. In addition, we assume little to no remanent magnetization and little to no self-demagnetisation. Important remanence and selfdemagnetization can be observed in the vicinity of magnetite-rich IOCG deposits (e.g., Anderson and Logan, 1992;Austin et al., 2013), but we consider that there is no indication of such features along L26. Further to this, we make this assumption for the sake of simplicity as the main object of this paper is the introduction of a new 440 sequential inversion workflow and to show that it is applicable to field data.
Under these premises, the features the magnetic data presents can be exploited to improve the image of the coverbasement interface when integrated with prior information about the thickness of cover. In this context, the role of magnetic data inversion constrained by MT is therefore multiple:  to constrain the depth and extent of the magnetic anomalies and refine their geometry; 445  to analyse the compatibility between the constraints derived from MT and the magnetic data and resolve some small-scale structures not defined by the MT constrains;  to reduce the interpretation uncertainty of the cover-basement interface;  to propose new scenarios in relation to the composition of the basement (in terms of its magnetic susceptibility) and structure (through its lateral variations). 450 The depth of the cover-basement interface probability shown in Figure 10c is used to derive the domains required by the spatially varying bound constraints using in magnetic inversion.

Magnetic data preparation and extraction of prior information
We use the gridded reduced-to-pole (RTP) magnetic data from the Geological Survey of Queensland shown in 455 Figure  10 (https://geoscience.data.qld.gov.au/dataset/ds000018/resource/91106497-d463-4b83-8b01-1c5539ab40b1, last accessed on 27/07/2021). Prior to the 2D inversion of the data along the line L26, we manipulate and reformat the data. To account for variations in the measurements in the vicinity of the line, we extract data from a 800 meters wide band around the profile (L26) (Figure 10a), as shown in more details Figure   10a. To obtain data corresponding to a 2D rectilinear profile, we then calculate the weighted average of this subset 460 of the dataset by assigning weights inversely proportional to the square of the distance of the measurement to L26, as illustrated in Figure 10b.

Inversion setup and results
To reduce computing time, we truncate the sensitivity matrix of the magnetic data inverse problem under the assumption that model-cells beyond a given distance exert a negligible influence on the forward calculated data.
The sensitivity domain is reduced to a 25 km-radius cylinder of infinite height and depth around each data point in a moving sensitivity domain approach (see examples in Čuma et al., 2012;Čuma and Zhdanov, 2014) that 485 assumes negligible contribution of the models-cells beyond a certain distance from the measurement's location.
Using a similar approach, the results of Wilson et al. (2011), indicate that using a 25 km radius may result in approximately 98% accuracy in the calculation of the response, which we assume to suffice for the purpose of our application example. It our case, it allows a reduction of 67% of the size of the sensitivity matrix. We performed inversion using 12 threads on an Intel(R) Xeon(R) CPU E5-2630. 490 To examine the impact of different type of constraints, we first perform inversions using minimum prior information and successively increase the amount of prior information from unconstrained inversion to using MTderived intervals for multiple bound constraints. The inversions we run consist of the following cases: (1) constrained by homogenous smoothness constraints; (2) constrained by homogenous smoothness constraints with lower and upper bound constraints; 495 (3) constrained by homogenous smoothness constraints with global, multiple bound constraints; (4) constrained by homogenous smoothness constraints with local, multiple bound constraints derived from MT.
The results for inversion cases (1) through (4) are shown in Figure 12 The inversions reached a satisfactory data fit, exception made for the constrained inversion 4 (see the data fit in 505 Figure 12e). In that case a significant underfit of the magnetic data is observed within certain areas, which point to an incompatibility between the magnetic data and the constraints applied. Four areas in the central part of the model are slightly underfit, on length scales < 300 meters. On the eastern part of the profile, from 479 km East to the most eastern part of the profile, an important underfit is observed. At this stage, this data misfit can indicate that the constraints used are not appropriate. This could be due to an inexact positioning at depth of the structural 510 constrains, or to a change in the petrophysical behaviour of the basement in certain areas, which would link differently the electrical properties of the depth to basement constrains to their magnetic properties. We propose a 5 th inversion case where we adjust the bounds manually to examine hypotheses relaxing the constraints derived by the combination of MT inversions and the magnetic susceptibility of rocks in the area.
From Figure 12, we identify five main areas where hypotheses made for the utilisation of MT-derived domains 515 need to be adjusted. In each case, the domain allowing sedimentary units may be deeper than expected or the basement may be less susceptible. We test the plausibility of such alternative scenarios and adapt the MT-derived domains by adjusting the domains. We increased the depth of the non-sedimentary (i.e., basement) units in the https://doi.org/10.5194/se-2021-124 Preprint. Discussion started: 3 November 2021 c Author(s) 2021. CC BY 4.0 License. eastern part of the model and between the areas marked by dashed lines in Figure 12, where the basementsedimentary interface may be deeper than anticipated or the basement characterized by less susceptible rocks. The 520 domains we use after adjustment are shown in Figure 13a, and inversion results in Figure 13b and c, respectively. Beyond the possibility to revisit hypotheses made at earlier stages of the workflow, we get insights into the structure and magnetic susceptibility of the basement. While electrical conductivity and magnetic susceptibility may be sensitive to change in rock type, there are scenarios where they exhibit differing sensitivity to texture and grain properties, respectively. For instance, metamorphism and alteration might affect electrical conductivity and 535 magnetic susceptibility differently (Clark, 2014;Dentith et al., 2020). Under these circumstances, our results can provide insights into plausible geological processes given sufficient prior geological information about the deformation history.

Interpretation
From a multi-physics modelling point of view, the results shown in the previous Sect. shows a general agreement 540 between the MT-derived constraints and the magnetic data. However, the results also point to incompatibilities in a few parts of the model. We identified two major areas where incompatibility occurs: 1) a smaller inconsistent areas in the western part of the survey; 2) a large inconsistent area east of the Mount Margaret Fault.
We interpret these incongruencies as being mainly due to the different sensitivities of the two geophysical methods 545 to different geological features and to the petrophysical variability of the basement in the area.
The greater depth extent of some of the lower magnetic susceptibility zones required by the magnetic data in the western part of the survey suggests that the depth to magnetic source is greater than suggested by the constraints.
Adjustments to the constraints allowed a better data fit. A low magnetic response between kms 460 and 470 East ( Figure 10) is assumed to be the response to low magnetic susceptibilities and is interpreted to be granitic 550 intrusions of the William-Naraku Supersuite (Dhnaram and Greenwood, 2013). The presence of such intrusions offers a plausible explanation for the discrepancies between the magnetic and MT modelling. On the one hand, the MT data modelling might not able to distinguish between an electrically resistive basement and an electrically resistive intrusion, while the magnetic data modelling could not distinguish between a non-magnetic cover and a non-magnetic intrusion. On the other hand, magnetic data inversion can differentiate the low susceptibility 555 intrusion from the higher susceptibility volcanic rocks, and the MT data is sensitive to the basal cover interface above both the volcanic rock and the intrusion. The constrained inversion permits detection of the lateral extent of the intrusion while estimating cover thickness. While detailed modelling of higher resolution data would be required to refine the geometry of these intrusive bodies, our modelling suggests that the intrusion could be modelled as several smaller intrusions. 560 East of the Mount Margaret Fault, the incompatibility between the original MT-derived constraints and the magnetic data points to regional scale structures. Drill hole observations indicate basement not exceeded 350 meters depth. If we assume a high-susceptibility basement, which is common to the whole area (Dhnaram and Greenwood, 2013), the magnetic model requires a very thick non-magnetic cover layer to reconcile the data which is incompatible with our geological knowledge of the area. In that case, we need to reconsider our definition of 565 the basement. The north-trending Mount Margaret Fault (see Figure 9)

Discussion
We have presented a workflow for sequential joint modelling of geophysical data, and applied it to synthetic and field measurements. In this study we used constraints in the form of interface probabilities derived from a probabilistic workflow driven by MT data, but it is general in nature and is not limited to a particular geological 575 or geophysical modelling method to generate the inputs. This has allowed us to report the utilisation of the ADMM algorithm to constrain magnetic data inversion for the first time This workflow presents several advantages. It is computationally inexpensive by use of standalone inversions.
The inversion of the MT dataset used to derive the constraints was performed only once. Then, a series of constrained magnetic inversions was run to test different geophysical and petrophysical hypotheses. It shows the 580 example of a fast and flexible approach to test different structural and petrophysical assumptions while modelling data sensitive to different physical parameters. It allows to focus the modelling efforts on survey-specific features (anomalies, geological structures) when appropriate petrophysical information is available. However, as with generalizable methods, strengths become limitations under certain circumstances. For instance, in the case of MT and magnetic data inversions as proposed in this work, the electrical resistivity and magnetic susceptibility for the 585 rock types of interest is dependent on a range of factors and processes (such as porosity, permeability, rock alteration, etc.) such that their correlation may be case-dependant (see Dentith et al., 2020;Dentith and Mudge, 2014). While we may surmise that it remains reasonable to assume the existence of such correlation in hard rock scenarios, it may not always hold in basin environments. For example, one can easily think of a basin exploration case where electrical resistivity increases rapidly with increasing hydrocarbon concentration in reservoirs, while 590 changes in magnetic susceptibility might make the use of magnetic data inversion redundant. In such case, other property pairings can be considered, such as electrical resistivity and seismic attributes (see examples of Le et al., 2016;Tveit et al., 2020, who use seismic inversion to extract prior information for CSEM inversion). Further to this, the utilisation of magnetic data inversion for the deeper part of the crust is limited to depths shallower than the Curie point's (typically from approximately 10 to a few 10s of kilometres under continents). For deeper 595 imaging of the crust, the workflow we propose may be suited to the utilisation of gravity data with MT.
An assumption worth examination is whether the study area is adequately represented by two geological domains.
In the cases we investigated, these domains were defined by the probability of observing only two rock classes (basement and non-basement). While this assumption reduces the risk of misinterpretation as no hypotheses are made to distinguish between the different sedimentary units or rocks of different nature in the basement, it also 600 then limits the interpretations that can be made from the results presented. We expect that provided that the rock units present 'desirable' features, i.e., distinctive magnetic susceptibility and resistivities (or other properties depending on the geophysical techniques considered) several rock types can be considered in the modelling. Such discriminative aspects of the petrophysics needs to be ascertained while defining the number of distinctive domains that may be present in the study area. This necessarily requires robust petrophysical data to be available 605 given the strong constraint that these domains impart on inversion. However, in the absence of petrophysical data or the number and character of geological domains, the magnitude of misfit can inform whether a proposed number of domains is plausible, driving data acquisition or refinement of the conceptual geological model.
The application case was performed in 2D to illustrate the workflow. Extending the presented work to large scale 610 problems in 3D is straightforward as the inversion methods employed in this study were designed for 3D modelling. The 1D MT modelling and interpolation schemes present excellent scalability. The Tomofast-x engine Ogarko et al., 2021b) is implemented using 3D grids. It presents good scalability and it offers the possibility to reduce the size of the computation domain to save memory when calculating the sensitivity matrix in the same fashion as Čuma et al. (2012), and Čuma and Zhdanov (2014), for large-scale potential field 615 data modelling. On this basis, the workflow presented here might be applicable to large datasets using, for publicly available potential field data. Studies of this kind would span continents for crustal investigations to 620 discover tectonic discontinuities and terrane interfaces. Current investigations on Tomofast-x comprise the application of wavelet compression operators to accelerate the inversion in the same way as Li and Oldenburg (2003) and Martin et al. (2013) and the development of joint inversions using the ADMM constraints for multiple bound constraints.
Another straightforward extension of the workflow is the use of gravity data simultaneously with, or instead of, 625 magnetic data since it is already implemented in Tomofast-x . Giraud et al. (2020) presented a synthetic MT-constrained gravity inversion, using a similar workflow as the one presented here. This would be of particular interest in the Cloncurry region (Queensland, Australia), where for instance, Moorkamp (2021) recently investigated the joint inversion of gravity and MT data, and where our workflow could be applied using the MT modelling results of Seillé et al. (2021). 630 From a geophysical point of view, magnetic inversion is affected by the non-uniqueness of the solution to the inverse potential field problem despite prior information and constraints being used. The workflow could be improved by using a series of models representative of the geological archetypes that can be derived from the ensembles of 1D MT models. Geological archetypes are distinctly different structural configurations (or topologies) that plausibly exist for a given location with available data (Pakyuz-Charrier et al., 2019, Wellmann 635 andCaumon, 2018). Identification of the archetypes could be achieved from the ensemble of geological model realisations in the same spirit as Pakyuz-Charrier et al. (2019), who use a Monte Carlo approach to generate a range of topologies which are then examined for distinct clusters which represent the archetypes.
From a methodological point of view, it could be argued that simultaneous joint geophysical inversion combining structural and petrophysical constraints might outperform the workflow we proposed here. However, this would 640 make the modelling process more demanding combined with limitations based on cases where determining the causitive relationships between petrophysics supporting joint approaches poses a challenge. The workflow we propose here presents a few advantages over a joint inversion scheme, in the sense that it does not require both datasets to be inverted simultaneously under a defined set of petrophysical and/or structural constraints. The time required to run a joint inversion being limited by the running time of the more computationally expensive 645 technique, it can limit the range of tests to be performed. In this study, we could run rapidly many 2D constrained magnetic inversions, even if the 1D probabilistic inversions of the MT data (and posterior fusion) required significant longer running time compared to the 2D constrained magnetic inversion. This point would particularly https://doi.org/10.5194/se-2021-124 Preprint. Discussion started: 3 November 2021 c Author(s) 2021. CC BY 4.0 License. be relevant in the case of large 3D datasets. This approach may represent a step in the modelling workflow which that is useful to explore, understand and refine structural and petrophysical relationships between different 650 physical parameters before undertaking more demanding joint inversions.
In the field application case presented here, the probabilistic depth to basement map to derive constraint assumed lateral continuity of the depth to basement estimates at a large scale, not accounting for small-scale lateral variations. Thus, uncertainty for depth to basement may be underestimated at some locations, in particular in between MT sites as shallow depths. In such cases, the existence of incompatibilities between MT-derived 655 constraints and the magnetic data might require reconsidering the spatial continuity assumptions taken during the calculation of the probabilistic depth to basement surface. Extensions of this work may be devised to alleviate some of the limitations of the workflow. For instance, magnetic susceptibility from the inversion of magnetic data could be mapped back to a resistivity model in order to calculate forward MT data for validation (dashed line in Figure 1), or to constrain the next cycle of MT inversions in the case the workflow is extended to cooperative joint 660 inversion. It would also be straightforward to use to a level-set inversion that can consider an arbitrary number of geological units (e.g., Giraud et al., 2021) using MT modelling as a source of prior information and constraints.
We have used hard bounds using the ADMM algorithm, which can easily be complemented or replaced by the use of multi-modal petrophysical distributions as available in Tomofast-x (e.g., mixture models as in Giraud et al., 2017Giraud et al., , 2019 as an alternative. 665

Conclusion
We have introduced, tested on synthetic, and applied to field data a cooperative inversion scheme for the integration of MT and magnetic inversions. We have shown that despite its simplicity, the workflow we propose efficiently leverages the complementarities between the two methods and has the capability to improve our understanding of the cover-basement interface and of the basement itself. We have tested our workflow on a 670 synthetic study that illustrates the flexibility of the method and the different possibilities our workflow offers as well as their limitations. In the field application case (Cloncurry area, Queensland), we have shown how the quantitative integration of MT and magnetic data may bring insightful results on geological structural and petrophysical aspects, opening up new avenues for interpretations of the geology of the area and prompting future works. 675 7 Authors' contribution JG and HS designed the study with inputs from the rest of the authors. JG and HS and generated the magnetic and MT synthetic data, respectively. JG performed the magnetic data inversions. JG redacted the manuscript with inputs from HS who is the main contributor to MT-related aspects, and comments from the rest of the authors. HS The updated model +1 is then calculated by solving the inverse problem using the LSQR algorithm of Paige 875 and Saunders (1982). We refer the reader to Ogarko et al. (2021a) for more details. We illustrate the application of such projection using two intervals in Figure A 1.