## 2 Methodology
### Bayesian Inference and the Kalman Ensemble Generator  
In Bayesian inversion approaches, all of the available information is translated into a probability distribution. Bayesian updates modify prior information by relevant measurement data producing posterior probability distributions for model parameters. The prior information is given as a probability density function (PDF) \(\rho(\mathbf{m})\) for the model parameters \(\mathbf{m}\in\mathbb{R}^{n_{par}}\). Inserting \(\mathbf{m}\) into the relevant geophysical forward model \(g\), and combining with the observations \(\mathbf{d}\in\mathbb{R}^{n_{obs}}\), we derive the likelihood function \(\rho(\mathbf{d}|\mathbf{m})\), the probability for a set of observations given the parameters. Here, \(n_{par}\) and \(n_{obs}\) denote the number of parameters and number of observations, respectively. [PERSON] theorem defines the posterior probability density function in terms of the PDFs discussed above:  
\[\rho(\mathbf{m}|\mathbf{d})=\kappa\cdot\rho(\mathbf{d}|\mathbf{m})\rho( \mathbf{m}), \tag{1}\]  
where random variables left of a vertical bar are conditional to a set of values for the random variables on the right side of the bar and \(\kappa\) denotes a normalization factor. In joint inversions, \(\mathbf{d}\) consists of data of two different types of geophysical measurements. Assuming that the measurement errors between the two data sets are uncorrelated, [PERSON] et al. [23] point out that the likelihood in a joint Bayesian framework can simply be given as the product of the individual likelihoods:  
\[\rho(\mathbf{d}|\mathbf{m})=\rho(\mathbf{d}_{\text{res}}|\mathbf{m})\rho( \mathbf{d}_{\text{$fdem$}}|\mathbf{m}) \tag{2}\]  
with a joint data vector  
\[\mathbf{d}= \begin{bmatrix}\mathbf{d}_{\text{res}}\\ \mathbf{d}_{\text{$fdem$}}\end{bmatrix}  \tag{3}\]  
for VES data \(\mathbf{d}_{\text{res}}\) and FDEM data \(\mathbf{d}_{\text{$fdem$}}\). Dealing with high-dimensional non-linear problems, no analytical solution is available to solve Equation (1). Instead, the posterior needs to be numerically sampled. Most commonly, this is done using Markov chain-Monte Carlo (MCMC) methods applying a Metropolis- or Gibbs-sampler. Using an infinite number of samples, this approach converges to the exact solution of the Bayesian update problem [24]. However, in practice, the computational resources are limited and the applicability of MCMC methods is therefore restricted. Computational restrictions are particularly pressing when expensive forward models are necessary and large data sets are used, which is usually the case in geophysical applications. Often, an MCMC sampler is no longer feasible for computing solutions to such problems as no convergence of the chain can be achieved within the number of affordable forward model runs. Here, we will apply a specific Monte Carlo implementation of the Bayesian update problem: the ensemble Kalman filter (EnKF, [15]). The EnKF was introduced as an efficient alternative to the original Kalman filter [25] for solving large-dimensional and non-linear problems in data assimilation [26]. In the EnKF context, the set of Monte Carlo samples is called an ensemble. The ensemble formulation makes the EnKF more efficient than the Kalman filter for two reasons: (1) in contrast to the the Kalman filter, the covariance equation can be replace by sample covariance and (2) Jacobian computations are unnecessary avoiding many calls of the forward model. Instead of the Jacobian, the full forward response \(g\) is computed for all ensemble realizations:  
\[\mathbf{d}^{i}=g(\mathbf{m}_{prior}^{i})+\epsilon_{err},\qquad i\in\{1,...,n_{ ens}\}, \tag{4}\]  
where \(n_{ens}\) is the size of the ensemble, \(\epsilon_{err}\) is the measurement error, and \(m_{i}\) and \(d_{i}\) are realisations of the Gaussian prior PDF and Gaussian data PDF, respectively. The ensemble is generated by drawing multivariate normal random numbers from the multivariate Gaussian prior using the Matlab function'mvnrnd' [27].