# Robust 3D Joint Inversion of Gravity and Magnetic Data: A High-Performance Computing Approach
## 2 Methods
### Regularized Least Squares 3D Joint Inversion with Gramian-Based Constraint
The linear operators L, M, and N correspond to the Laplacian \(L=\left(\frac{\partial^{2}f}{\partial x^{2}}+\frac{\partial^{2}f}{\partial y^ {2}}+\frac{\partial^{2}f}{\partial z^{2}}\right)\) which smoothes the distribution of properties in the solutions, \(M=\left(\frac{\partial f}{\partial x}cos\phi\text{ cos}\theta+\frac{\partial f}{ \partial y}cos\phi\text{ sin}\theta+\frac{\partial f}{\partial z}sin\phi\right)\), which is the directional derivative of the strike \(\theta\) and pitch \(\phi\) angles and whose function is to specify a preferential structural direction, as well as the vertical derivative \(N=\left(\frac{\partial f}{\partial z}\right)\), which depending of its regularization parameter sign favors vertical or horizontal properties distributions in the subsurface. These matrices were computed from the linearization of the 3D domain and discrete numerical derivatives using central differences and forward and backward differences depending on which position of each node in the 3D domain (Figure 1). This allowed us to obtain the curvature and derivatives of the domain without the need to consider extra cells or reduce the domain. The operators \(A^{(1)}\) and \(A^{(2)}\) were calculated using analytical solutions for straight prisms following the development of [77,78] for gravimetry and magnetometry, respectively (See Appendix C). The Gramian \(G\)[73] is defined as the determinant of the Gram matrix for the elements \(m^{(1)}\) and \(m^{(2)}\). The equivalence of the Gramian with respect to the L2 norm of crossed gradients is detailed in Appendix B of [76]. \[G\Big{(}\
abla m^{(1)},\
abla m^{(2)}\Big{)}=det\left| \begin{pmatrix}\
abla m^{(1)},\
abla m^{(1)}\\ \
abla m^{(2)},\
abla m^{(1)}\end{pmatrix} \right.\quad \begin{pmatrix}\
abla m^{( 1)},\
abla m^{(2)}\\ \
abla m^{(2)},\
abla m^{(2)}\end{pmatrix} \right|=\|\
abla m^{(1)}\times \
abla m^{(2)}\|^{2} \tag{4}\]  
\[G\Big{(}\
abla m^{(2)},\
abla m^{(1)}\Big{)}=det\left| \begin{pmatrix}\
abla m^ {(2)},\
abla m^{(2)}\\ \
abla m^{(1)},\
abla m^{(2)}\end{pmatrix} \right.\quad \begin{pmatrix}\
abla m^ {(2)},\
abla m^{(1)}\\ \
abla m^{(1)},\
abla m^{(1)}\end{pmatrix} \right|=\|\
abla m^{(2)}\times \
abla m^{(1)}\|^{2} \tag{5}\]  
From Equation (2), we obtain [PERSON]'s equation, which expresses the parameters \(m\) at which the parametric functional reaches a minimum (See Appendix B):  
\[ \begin{array}{c}I(m)=A^{T}W_{d}^{2}A(m)-A^{T}W_{d}^{2}d+\alpha L^{T}L(m)+ \beta W_{m}^{2}(m)-\beta W_{m}^{2}(m_{apr})\\ +\gamma M^{T}M(m)+\eta N^{T}N(m)+\mu I_{G}\end{array}  \tag{6}\]  
Regrouping the terms of Equation (6) as follows will be useful for their computational implementation:  
\[I(m)=(I_{A})m-I_{B}+\mu I_{G} \tag{7}\]  
where:  
\[I_{A}=A^{T}W_{d}^{2}A+\alpha L^{T}L+\beta W_{m}^{2}+\gamma M^{T}M+\eta N^{T}N \tag{8}\]  
\[I_{B}=A^{T}W_{d}^{2}d+\beta W_{m}^{2}(m_{apr}) \tag{9}\]  
\[I_{G}^{(1)}=\
abla\,\cdot\Big{(}\Big{\langle}\,\
abla m_{n}^{(2)},\
abla m_{n} ^{(1)}\Big{\rangle}\
abla m_{n}^{(2)}\Big{)}-\
abla\cdot(\|\
abla m_{n}^{(2)} \|^{2}\
abla m_{n}^{(1)}\Big{)} \tag{10}\]  
\[I_{G}^{(2)}=\
abla\,\cdot\Big{(}\Big{\langle}\,\
abla m_{n}^{(1)},\
abla m_{ n}^{(2)}\Big{\rangle}\
abla m_{n}^{(1)}\Big{)}-\
abla\cdot(\|\
abla m_{n}^{(1)} \|^{2}\
abla m_{n}^{(2)}\Big{)} \tag{11}\]  
\[I_{G}^{(1)}=\
abla\cdot\Big{(}\Big{\langle}\,\
abla m_{n}^{(2)},\
abla I_{n} ^{(1)}\Big{\rangle}\
abla m_{n}^{(2)}\Big{)}-\
abla\cdot(\|\
abla m_{n}^{(2)} \|^{2}\
abla I_{n}^{(1)}\Big{)} \tag{12}\]  
\[J_{G}^{(2)}=\
abla\cdot\Big{(}\Big{\langle}\,\
abla m_{n}^{(1)},\
abla I_{n} ^{(2)}\Big{\rangle}\
abla m_{n}^{(1)}\Big{)}-\
abla\cdot(\|\
abla m_{n}^{(1)} \|^{2}\
abla I_{n}^{(2)}\Big{)} \tag{13}\]  
\[I_{G}= \begin{bmatrix}I_{G}^{(1)}\\ I_{G}^{(2)}\end{bmatrix} \ \ ;\ \ J_{G}= \begin{bmatrix}I_{G}^{(1)}\\ I_{G}^{(2)}\end{bmatrix}  \tag{14}\]  
Details of the expressions \(I_{G}\,y\,\,\,J_{G}\) are given in Appendix C.  
Figure 1: 3D mesh illustrating the cells involved in the computation of discrete numerical derivatives depending on the position of the central node.